text:
"""
# Method p62

Plane strain bearing capacity analysis of an elastic-plastic
(von Mises) material using 8-node rectangular quadrilaterals.

Viscoplastic strain method.

No global stiffness matrix assembly.

Diagonally preconditioned conjugate gradient solver.

### Constructors
```julia
p62(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}` : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Structural element
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Floalt64}                       : x-coordinate vector
* y_coords::FloatRange{Floalt64}                       : y-coordinate vector
* thickness:: Float64                                  : Thickness of plate
* tol::Float64                                         : Convergence tolerance
* qincs::Vector{Float64}                               : Incremental load steps
```

### Optional additional data dictionary keys
```julia
* limit = 250                  : Iteration limit
* penalty = 1e20               : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}         : Element material vector if np_types > 1
```

### Return values
```julia
* (g_coord, g_num, disp)        : where:
                                    g_coord  : Coordinates
                                    g_num    : Node numbering
                                    disp     : Matrix of displacements
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Plane                         : Help on a Plane structural element
?FiniteElement                 : List finite element types
?Quadrilateral                 : Help on Quadrilateral finite element
```
"""
function p62(data::Dict)

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
    @assert typeof(struc_el) <: StructuralElement
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  if struc_el.axisymmetric
    nst = 4
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  tol = 1.0e-10
  if :tol in keys(data)
    tol = copy(data[:tol])
  end

  limit = 100
  if :limit in keys(data)
    limit = copy(data[:limit])
  end

  if :cg_tol in keys(data)
    cg_tol = data[:cg_tol]
  else
    cg_tol = 1.0e-5
  end

  if :cg_limit in keys(data)
    cg_limit = data[:cg_limit]
  else
    cg_limit = Int(1e3)
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  g_coord = zeros(ndim,nn)
  if :g_coord in keys(data)
    g_coord = data[:g_coord]'
  end

  g_num = zeros(Int, fin_el.nod, nels)
  if :g_num in keys(data)
    g_num = reshape(data[:g_num]', fin_el.nod, nels)
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  eps = zeros(nst)
  stress = zeros(nst)
  bload =zeros(ndof)
  eload =zeros(ndof)
  erate = zeros(nst)
  evp = zeros(nst)
  devp = zeros(nst)
  m1 = zeros(nst, nst)
  m2 = zeros(nst, nst)
  m3 = zeros(nst, nst)
  flow = zeros(nst, nst)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)

  @assert :qincs in keys(data)
  qincs = data[:qincs]

  println("There are $(neq) equations.")

  loads = zeros(Float64, neq+1)
  bdylds = zeros(Float64, neq+1)
  oldis = zeros(Float64, neq+1)
  totd = zeros(Float64, neq+1)
  p = zeros(nst, struc_el.nip, nels)
  x = zeros(neq+1)
  u = zeros(neq+1)
  xnew = zeros(neq+1)
  diag_precon = zeros(neq+1)
  d = zeros(neq+1)

  # Find global array sizes
  for iel in 1:nels
    geom_rect!(fin_el, iel, x_coords, y_coords, coord, num, struc_el.direction)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
  end

  tensor = zeros(nst, struc_el.nip, nels)
  storkm = zeros(ndof, ndof, nels)
  #tload = zeros(neq+1)
  #dtemp = zeros(nn)
  #dtel = zeros(fin_el.nod)
  #epsi = zeros(nst)
  dt = 1.0e15
  ddt = 0.0
  sigm = 0.0
  dsbar = 0.0
  lode_theta = 0.0

  sample!(fin_el, points, weights)

  # element stiffnes integration, storage and conditioner

  for iel in 1:nels
    ddt = 4.0*(1.0+prop[etype[iel], 3])/(3.0*prop[etype[iel], 2])
    if ddt < dt
      dt = ddt
    end
    deemat!(dee, prop[etype[iel], 2], prop[etype[iel], 3])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      beemat!(bee, deriv)
      km += (bee')*dee*bee*detm*weights[i]
    end
    storkm[:,:,iel] = km
    for k in 1:ndof
      diag_precon[g[k]+1] = diag_precon[g[k]+1] + km[k,k]
    end
  end

  diag_precon = 1 ./ diag_precon

  # Read load weightings

  loaded_nodes = 0
  node = Int[]
  val = Array{Float64, 2}
  if :loaded_nodes in keys(data)
    loaded_nodes = size(data[:loaded_nodes], 1)
    node = zeros(Int, loaded_nodes)
    val = zeros(loaded_nodes, size(data[:loaded_nodes][1][2], 2))
  end

  if :loaded_nodes in keys(data)
    for i in 1:loaded_nodes
      node[i] = data[:loaded_nodes][i][1]
      val[i,:] = data[:loaded_nodes][i][2]
    end
  end

  nf1 = deepcopy(nf) .+ 1

  # Load increment loop

  #println("\nstep     load        disp       iters    cg/plastic iters")

  converged = false
  ptot = 0.0
  iters = 0
  cg_tot = 0
  diag_precon[1] = 0.0
  iy = 0

  # DataFrame arrays

  loaddt = Float64[]
  dispdt = Float64[]
  itersdt = Int[]
  ratiodt = Float64[]

  #for iy in 1:1
  for iy in 1:size(qincs, 1)
    ptot = ptot + qincs[iy]
    iters = 0
    bdylds = zeros(Float64, neq+1)
    evpt = zeros(nst, struc_el.nip, nels)
    cg_tot = 0
    cg_converged = false

    # Plastic iteration loop

    while true
      iters += 1
      loads = zeros(Float64, neq+1)
      for i in 1:loaded_nodes
        loads[nf1[:,node[i]]] = val[i,:] * qincs[iy]
      end
      loads += bdylds
      d = diag_precon .* loads
      p = d
      x = zeros(neq+1)
      cg_iters = 0

      # Pcg equation solution
      tmpconvcrit = 0.0

      while true
        cg_iters += 1
        u = zeros(neq+1)

        for iel in 1:nels
          g = g_g[:, iel]
          km = storkm[:, :, iel]
          u[g .+ 1] += km * p[g .+ 1]
        end
        up = dot(loads, d)
        alpha = up ./ dot(p, u)
        xnew = x + p .* alpha
        loads -= u .* alpha
        d = diag_precon .* loads
        beta = dot(loads, d) ./ up
        p = d + p * beta
        #tmpconvcrit =  maximum(abs.(xnew-x))/maximum(abs.(xnew))
        cg_converged = checon(xnew, x, cg_tol)
        #iy < 2 && println([iters cg_iters cg_tot+cg_iters tmpconvcrit])
        if cg_converged || (cg_iters == cg_limit)
          break
        end
      end
      cg_tot += cg_iters
      loads = xnew
      loads[1] = 0.0

      # Check plastic convergence

      tmp = maximum(abs.(loads-oldis))/maximum(abs.(loads))
      converged = checon(loads, oldis, tol)
      if iters == 1
        converged = false
      end
      if converged || (iters == limit)
        bdylds = zeros(Float64, neq+1)
      end

      # Go round the Gauss point

      for iel in 1:nels
        deemat!(dee, prop[etype[iel], 2], prop[etype[iel], 3])
        num = g_num[:, iel]
        coord = g_coord[:, num]'
        g = g_g[:, iel]
        eld = loads[g .+ 1]
        bload = zeros(ndof)
        for i in 1:struc_el.nip
          shape_der!(der, points, i)
          jac = der*coord
          detm = det(jac)
          jac = inv(jac)
          deriv = jac*der
          beemat!(bee, deriv)
          eps = bee*eld
          eps -= evpt[:, i, iel]
          sigma = dee*eps
          stress = sigma + tensor[:, i, iel]
          (sigm, dsbar, lode_theta) = invar(stress, sigm, dsbar, lode_theta)

          # Check whether yield is violated

          f = dsbar - sqrt(3.0)*prop[etype[iel], 1]
          if converged || iters == limit
            devp = deepcopy(stress)
          else
            if f >= 0.0
              dq1 = 0.0
              dq2 = 3.0/2.0/dsbar
              dq3 = 0.0
              (m1, m2, m3) = formm!(stress, m1, m2, m3)
              flow = f*(m1*dq1 + m2*dq2 + m3*dq3)
              erate = flow*stress
              evp = erate*dt
              evpt[:,i,iel] += evp
              devp = dee*evp
            end
          end
          if f >= 0.0 || converged || iters == limit
            eload = devp' * bee
            bload += eload' .* detm .* weights[i]
          end

          # Update the Gauss point stresses

          if converged || iters == limit
            tensor[:,i,iel] = stress
          end
        end

        # Compute the total bodyloads vector

        bdylds[g .+ 1] += bload
        bdylds[1] = 0.0
      end
      if converged || iters == limit
        break
      end
    end
    totd += loads
    totdstr = @sprintf("%+.4e", totd[nf1[2, node[1]]])

    append!(loaddt, [ptot])
    append!(dispdt, [totd[nf1[2, node[1]]]])
    append!(itersdt, [iters])
    append!(ratiodt, [round.(cg_tot/iters, digits=2)])

    if iy < 10
      println(" $(iy)       $(ptot)    $(totdstr)    $(iters)       $(round.(cg_tot/iters, digits=2))")
    else
      println("$(iy)       $(ptot)    $(totdstr)    $(iters)       $(round.(cg_tot/iters, digits=2))")
    end

  end
  println()

  dt = DataFrame()
  dt[!, :loads] = loaddt
  dt[!, :disp] = dispdt
  dt[!, :iters] = itersdt
  dt[!, :ratio] = ratiodt
  dt
end

max_stars_repo_path: src/6 Material Nonlinearity/p62.jl

max_issues_repo_path: src/6 Material Nonlinearity/p62.jl

max_forks_repo_path: src/6 Material Nonlinearity/p62.jl


text:
module SteinDiscrepancy

    # We need to make sure the depsjl file is loaded so spanner .so is reachable
    depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    isfile(depsjl) ? include(depsjl) : error("SteinDiscrepancy not properly ",
        "installed. Please run\nPkg.build(\"SteinDiscrepancy\")")

    using MathProgBase.SolverInterface: AbstractMathProgSolver
    using JuMP
    using LinearAlgebra
    using Distributed

    import Base.rand
    import DataStructures
    import Distributions.TDist
    import Distributions.quantile

    export SteinDiscrete,
           # Main method for computing Stein discrepancy
           stein_discrepancy,
           # The stein_discrepancy aliases
           ksd,
           gsd,
           # Diffusion specific stein discrepancies
           langevin_graph_discrepancy,
           riemannian_langevin_graph_discrepancy,
           langevin_classical_discrepancy,
           langevin_kernel_discrepancy,
           # Result objects returned by affine stein discrepancy methods
           AffineDiscrepancyResult,
           # Result objects returned by langevin_kernel_discrepancy w/o checkpoints
           LangevinKernelResult,
           # Result objects returned by langevin_kernel_discrepancy w/ checkpoints
           LangevinKernelCheckpointsResult,
           # Function for computing spanner edges
           getspanneregdes,
           # Function computing univariate Wasserstein distancen
           wasserstein1d,
           # Function computing Wasserstein distance for discrete distributions
           wassersteindiscrete,
           # Function approximating Wasserstein distance for multivariate Wasserstein
           approxwasserstein,

           ### kernel exports
           SteinKernel,
           SteinParzenARKernel,
           SteinChampionLenardMillsKernel,
           SteinGaussianWeightedKernel,
           SteinGaussianPowerKernel,
           SteinGaussianUnboundedPowerKernel,
           SteinGaussianKernel,
           SteinMaternTensorizedKernel,
           SteinMaternRadialKernel,
           SteinPolyharmonicSplineKernel,
           SteinInverseMultiquadricKernel,
           SteinRectangularDomainMetaKernel,
           ### kernel utility functions
           k,
           gradxk,
           gradyk,
           gradxyk,
           k0

    ### distributions
    include("SteinDiscrete.jl")

    ### kernels
    include("kernels/SteinKernel.jl")
    include("kernels/SteinTensorizedKernel.jl")
    include("kernels/SteinGaussianWeightedKernel.jl")
    include("kernels/SteinGaussianPowerKernel.jl")
    include("kernels/SteinGaussianUnboundedPowerKernel.jl")
    include("kernels/SteinGaussianKernel.jl")
    include("kernels/SteinParzenARKernel.jl")
    include("kernels/SteinChampionLenardMillsKernel.jl")
    include("kernels/SteinPolyharmonicSplineKernel.jl")
    include("kernels/SteinMaternTensorizedKernel.jl")
    include("kernels/SteinMaternWeightedKernel.jl")
    include("kernels/SteinMaternRadialKernel.jl")
    include("kernels/SteinMaternPowerKernel.jl")
    include("kernels/SteinMaternUnboundedPowerKernel.jl")
    include("kernels/SteinInverseMultiquadricKernel.jl")
    # meta kernels
    include("kernels/SteinRectangularDomainMetaKernel.jl")

    ### discrepancy
    # Spanner files and utilities
    include("discrepancy/spanner.jl")
    include("discrepancy/l1spanner.jl")
    include("discrepancy/mincostflow.jl")
    # graph related discrepancies
    include("discrepancy/AffineDiscrepancyResult.jl")
    include("discrepancy/affine_graph_discrepancy.jl")
    include("discrepancy/langevin_graph_discrepancy.jl")
    include("discrepancy/riemannian_langevin_graph_discrepancy.jl")
    # classical related discrepancies
    include("discrepancy/affine_classical_discrepancy.jl")
    include("discrepancy/langevin_classical_discrepancy.jl")
    # kernel related discrepancies
    include("discrepancy/LangevinKernelBaseResult.jl")
    include("discrepancy/langevin_kernel_discrepancy.jl")
    # utilities / addons
    include("discrepancy/wasserstein.jl")
    # the main ONE
    include("discrepancy/stein_discrepancy.jl")
end

ext: jl


text:
#=
This example follows the approach by Luer et al. doi:10.1021/acs.jpclett.6b02704 to simulate
Frenkel excitons in TDBC J-aggregates.

Essentially it averages many different realisations of a J-aggregate, resulting in an
averaged, smooth lineshape.

=#

using MultidimensionalSpectroscopy, PyPlot, QuantumOptics, LinearAlgebra, FFTW, Colors,
        Printf, DelimitedFiles, Random, Combinatorics, Distributions

# make sure to set script directory as pwd()
cd(@__DIR__)

pygui(true)

calc_2d = false

try
    local m = @which dω(1)
    Base.delete_method(m)
catch
end

function dω(stdev)
    return rand_normal(0,stdev)
end

tlist = [0:5.:750;]
corr = zeros(length(tlist))

cmp = create_colormap("bright");

N = 70

#E = [2.41 * dE_n[i] for i in 1:N]
E_mon = 2.41
E_exc = 2.11
E_RWA = E_mon
E_RWA = E_exc
E     = E_mon - E_RWA

b_mon = NLevelBasis(2)
σ⁺ = transition(b_mon,2,1)
σ⁻ = transition(b_mon,1,2)

#rho0 = dm(tensor([nlevelstate(b_mon,1) for i in 1:N]))

B_exc = b_mon^N

H_n(n)  = embed(B_exc,n,E * σ⁺ * σ⁻) # Hamiltonian at site n => Hₙ

Σ⁺_n(n) = embed(B_exc,n,σ⁺)              # raising operator at site n
Σ⁻_n(n) = embed(B_exc,n,σ⁻)              # lowering operator at site n

#Σ⁺ = sum(Σ⁺_n.(1:N))
#Σ⁻ = sum(Σ⁻_n.(1:N))

N_deloc = 10
Δ = .299
#J = [-Δ / (2 * cos(pi /(N+1))) for i in 1:N]
J = -Δ / (2 * cos(pi /(N_deloc+1)))
J  = -0.15  # from Luer
J2 = -0.01 # from Luer
#J = 0
#J2 = 0
er = 3
e0 = 8.854e-21 #F/nm
ke = (4 * pi * er * e0)^-1
μd = 1.869e-20 # Cm

R = 1

#J = ke * (μd^2 / R^3)

#if N > 3
#    H = sum(H_n.(1:N)) +
#                sum(J[i] * (Σ⁺_n(i) * Σ⁻_n(i+1) + Σ⁻_n(i) * Σ⁺_n(i+1)) for i in 1:N-1)
#end

b = NLevelBasis(N)
p = transition(b,2,1)
H = dagger(p) * p * 0

rho0 = copy(H)
rho0.data[1,1] = 1

μ    = copy(H)

Ltemp   = copy(H)
L = []

for i in 2:N
    H.data[i,i] = E
    if i > 2
        H.data[i,i-1] = J
        H.data[i-1,i] = J
    end
    if i > 3
        H.data[i,i-2] = J2 # J2 coupling
        H.data[i-2,i] = J2
    end
    μ.data[1,i] = 1
    μ.data[i,1] = 1

    Ltemp.data[1,i] = 1
    append!(L,[copy(Ltemp)])
    Ltemp.data[1,i] = 0

end

#μ    = sum((Σ⁻_n(i) + Σ⁺_n(i)) for i in 1:N)

μ12  = rho0 * μ + μ * rho0
μ23  = copy(μ12)

rho1 = μ12 * rho0 * μ12
μ12  = μ12 / sqrt(tr(rho1))
rho1 = μ12 * rho0 * μ12
rho2 = μ23 * rho1 * μ23
#μ23  = μ23 / sqrt(tr(rho2))
rho2 = μ23 * rho1 * μ23

# TDBC details https://arxiv.org/pdf/1506.01321.pdf
γ01 = 0.00476 # eV ... from 1.15 # ps⁻¹ from photoluminescence
Γ01d = 0.017  # eV from 17 meV
Γ01 = γ01 / 2 + Γ01d
#L = Γ01 .* [Σ⁻_n(i) * (Σ⁺_n(i) * Σ⁻_n(i)) for i in 1:4]
L = Γ01 .* L

#=
idx = sortperm(real(diag((H).data)))
H.data = H.data[idx,idx]
Σ⁺.data = Σ⁺.data[idx,idx]
Σ⁻.data = Σ⁻.data[idx,idx]
rho0.data = rho0.data[idx,idx]
rho1.data = rho1.data[idx,idx]
rho2.data = rho2.data[idx,idx]
for i in 1:N
    #Σ⁺n[i].data = Σ⁺n[i].data[idx,idx]
    #Σ⁻n[i].data = Σ⁻n[i].data[idx,idx]
    L[i].data = L[i].data[idx,idx]
end

H_si, transf_mat_si, P_si, L_si, rho0_si, rho1_si, rho2_si, μ12_si, μ23_si, Σ⁺_si, Σ⁻_si =
    cmds.create_subspace([H],"si", L, rho0, rho1, rho2, μ12, μ23, Σ⁺, Σ⁻)

##

H_si = transf_mat_si * H_si[1] * transf_mat_si'
rho0_si = transf_mat_si * rho0_si * transf_mat_si'
rho1_si = transf_mat_si * rho1_si * transf_mat_si'
rho2_si = transf_mat_si * rho2_si * transf_mat_si'
μ12_si = transf_mat_si * μ12_si * transf_mat_si'
μ23_si = transf_mat_si * μ23_si * transf_mat_si'
Σ⁺_si = transf_mat_si * Σ⁺_si * transf_mat_si'
Σ⁻_si = transf_mat_si * Σ⁻_si * transf_mat_si'
=#

##
figure()

D = 1/100
tc = 100
gt = D^2 * tc^2 * (exp.(-tlist./tc) .+ tlist./tc .- 1)
#gt = tc^2 * (exp.(-tlist./tc) .+ tlist./tc .- 1)

F = L

M = 500
for i in 1:M
    println(i)

    global H_in_loop = copy(H)

    # intrachain disorder
    local dE_n = [dω(0.022) for i in 1:N]

    # interchain disorder
    local dE = dω(0.022)
    println(E + dE_n[1])

    #plot(collect(1:N),E .+ dE_n)

    for j in 2:N
        H_in_loop.data[j,j] = (H_in_loop.data[j,j] + dE_n[j-1] ) #+ dE
    end

    # introduce defect sites
    global N_defect = 0
    if rand() > 0.
        N_defect = 3
    else
        N_defect = 0
    end
    #global N_defect = rand(1:10)
    global N_defect = Int8(round(rand(Normal(5, .01))))
    if N_defect <= 0
        N_defect = 1
    end

    global ind_defect = rand(2:(N-1),N_defect)
    println(ind_defect)
    #if maximum(ind_defect) != N
        for j in ind_defect
            H_in_loop.data[j,j+1] = .6 * H_in_loop.data[j,j+1]
            H_in_loop.data[j+1,j] = .6 * H_in_loop.data[j+1,j]
            H_in_loop.data[j,j-1] = .6 * H_in_loop.data[j,j-1]
            H_in_loop.data[j-1,j] = .6 * H_in_loop.data[j-1,j]
            if j < N-1
                H.data[j,j+2] = .3 * H.data[j,j+2]
                H.data[j+2,j] = .3 * H.data[j+2,j]
            end
        end
    #else
    #end

    println("start corr")

    global corr = ( (timecorrelations.correlation(tlist, dense(rho0), H_in_loop, L, μ12, μ12)) + corr ) ./ 2

    #plot(tlist,corr)
    if i != M
        if i == 1
            plot([1, N], -1 .* [J, J],label="J")
        end
        scatter(collect(1:N) .+ rand(),dE_n .+ E_RWA,s=15,alpha=1/M*10,color="gray")
        scatter(ind_defect .+ .5,[.6 * -J for i in 1:length(ind_defect)], s=15,alpha=.01,color="red")
    elseif i == M
        scatter(collect(1:N) .+ rand(),dE_n .+ E_RWA,s=15,alpha=.25,color="black",label="site energies")
        scatter(ind_defect .+ .5,[.6 * -J for i in 1:length(ind_defect)], s=15,alpha=.25,color="red",label="defects")
    end

    if calc_2d
        if i < 3
            local T = [0]
            global spectra2d = Array{out2d}(undef, length(T))
            spectra2d[1] = make2Dspectra(tlist,[rho0, rho0],[H_in_loop, H_in_loop],[F, F],[μ12, μ12],[μ12, μ12],T[1],
                                                "lindblad";debug=false,use_sub=true,
                                                    t2coh=false,zp=10)
            spectra2d[1].ω = spectra2d[1].ω .+ E_RWA

            if i == 1
                global corr2d = spectra2d[1].corr
            else
                corr2d = (spectra2d[1].corr + corr2d) ./ 2
            end
        end
    end

end
ylim((0,3))
legend()

#corr /= M
corr =corr .* exp.(-gt)

freqRWA = exp.(-1im * E_RWA .* tlist) .* exp.(1im * E_RWA .* vcat(-reverse(tlist)[1:end-1],tlist))'

if calc_2d
    #corr2d = corr2d .* exp.(-gt) .* vcat(reverse(exp.(-gt))[1:end-1],exp.(-gt))'
    #spec2d  = corr2spec(corr2d .* freqRWA,10)

    figure()
    pcolor(real(spec2d))
end

zp = 11
corr    = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)

ω_exc, spec_exc = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true);

##

figure();
subplot(211)
plot(tnew,corr)
plot(tlist,exp.(-gt))

subplot(212)
plot(-1 .* (ω_exc.-E_RWA),spec_exc)

dataTDBC_mon = readdlm("TDBCmon_abs.dat", ',')
dataTDBC_aggr = readdlm("TDBCaggr_abs.dat", ',')

plot(dataTDBC_mon[:,1],dataTDBC_mon[:,2]/maximum(dataTDBC_mon[:,2]),"k:")
plot(dataTDBC_aggr[:,1],dataTDBC_aggr[:,2]/maximum(dataTDBC_aggr[:,2]),"r:")

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
using MatpowerCases, Base.Test

println("- extract_branch check")
for c in casenames()
    mpc = loadcase(c; describe=false)
    @assert isa(extract_branch(mpc), MatpowerCases.Branch)
    @assert isa(extract_branch(c), MatpowerCases.Branch)
end


text:
# And 2xN matrix of N landmark xy positions as variable nodes in factor graph
function addSimMapFG!(fg::FactorGraph, lms::Array{Float64,2})
    for i in 1:size(lms,2)
      newLandm!(fg,string('l',i), vectoarr2(lms[:,i]), 0.001*eye(2))
    end
    nothing
end

# assume all the landmarks are already loaded into the Ground Truth FG
function simOdo!(fgGT::FactorGraph, fg::FactorGraph, DX::Array{Float64,1};
    noiserate=2.0*[3e-2;3e-2;1.5e-3], driftrate=[0.0;0.0;0.0], detLM=Union{})
    prev, X, nextn = getLastPose2D(fg)
    addOdoFG!(fgGT, nextn, DX, 0.001*eye(3))

    r = norm(DX[1:2])
    xn = noiserate[1]*r
    yn = noiserate[2]*r
    thn = noiserate[3]*r
    cov = diagm([xn;yn;thn])
    DXn = DX + [xn*randn();yn*randn();thn*randn()] + r*driftrate
    addOdoFG!(fg, nextn, DXn, cov)

    return nextn
end

function truePredBR(fgGT::FactorGraph, fg::FactorGraph, ps::String, lm::String)
    trubr = predictBodyBR(fgGT, ps, lm)
    truA = [trubr[1]; trubr[2]]
    # Prediction of BR measurement
    prebr = predictBodyBR(fg, ps, lm)
    preA = [prebr[1];prebr[2]]
    return truA, preA
end

function showTruePredBR(fgGT::FactorGraph, fg::FactorGraph, ps::String, lm::String, cov::Array{Float64,2})
    truA, preA = truePredBR(fgGT, fg, ps, lm)
    measA = truA + [cov[1,1]*randn();cov[2,2]*randn()]
    mala = malahanobisBR(measA, preA, cov)
    println(ps,lm, ": true BR=$(round(truA,3)), pred BR=$(round(preA,3)), mala=$(round(mala,3))")
end

function crossMalaBR(fgGT::FactorGraph, fg::FactorGraph,
                      ps::String, lmT::String, lmE, cov::Array{Float64,2})
    truT, preT = truePredBR(fgGT, fg, ps, lmT)
    measT = truT + [cov[1,1]*randn();cov[2,2]*randn()]

    truE, preE = truePredBR(fgGT, fg, ps, lmE)
    malaX = malahanobisBR(measT, preE, cov)
    return malaX
end


text:
function Incrementor()
    return Incrementor(())
end

function Incrementor(arg0::jint)
    return Incrementor((jint,), arg0)
end

function Incrementor(arg0::jint, arg1::Incrementor_MaxCountExceededCallback)
    return Incrementor((jint, Incrementor_MaxCountExceededCallback), arg0, arg1)
end

function can_increment(obj::Incrementor)
    return jcall(obj, "canIncrement", jboolean, ())
end

function can_increment(obj::Incrementor, arg0::jint)
    return jcall(obj, "canIncrement", jboolean, (jint,), arg0)
end

function get_count(obj::Incrementor)
    return jcall(obj, "getCount", jint, ())
end

function get_maximal_count(obj::Incrementor)
    return jcall(obj, "getMaximalCount", jint, ())
end

function increment(obj::Incrementor)
    return jcall(obj, "increment", void, ())
end

function increment(obj::Incrementor, arg0::jint)
    return jcall(obj, "increment", void, (jint,), arg0)
end

function reset(obj::Incrementor)
    return jcall(obj, "reset", void, ())
end

function with_callback(obj::Incrementor, arg0::Incrementor_MaxCountExceededCallback)
    return jcall(obj, "withCallback", Incrementor, (Incrementor_MaxCountExceededCallback,), arg0)
end

function with_count(obj::Incrementor, arg0::jint)
    return jcall(obj, "withCount", Incrementor, (jint,), arg0)
end

function with_maximal_count(obj::Incrementor, arg0::jint)
    return jcall(obj, "withMaximalCount", Incrementor, (jint,), arg0)
end


text:
using BasisFunctions
using Overlap
using OffsetArrays

@testset "Overlap.jl" begin
    Rμ = SVector{3}(randn(3))
    Rν = SVector{3}(randn(3))
    ζμ = abs(randn())
    ζν = abs(randn())
    @test all(buildCartesianOverlaps(Rμ, ζμ, 0, Rν, ζν, 0)[0,0,:] .≈ 1.6392143992458958)
    @test begin
        S = view(buildCartesianOverlaps(Rμ, ζμ, 2, Rν, ζν, 3), :, :, 2)
        Target = OffsetArray([1.6392143992 -0.0126705658 0.9068452213 -0.0210272948;
                              0.1026258527  0.9059540195 0.0427569217  1.5035741369;
                              0.9131723512  0.1064784544 1.5065815043  0.1534566213], -1, -1)
        S ≈ Target
    end
end


text:
# Epi-composition, also known as infimal postcomposition or image function.

# This is the dual operation to precomposition, see Rockafellar and Wets,
# "Variational Analysis", Theorem 11.23.
#
# Given a function f and a linear operator L, their epi-composition is:
#
#   g(y) = (Lf)(y) = inf_x { f(x) : Lx = y }.
#
# Plugging g directly in the definition of prox, one has:
#
#   prox_{\gamma g}(z) = argmin_y { (Lf)(y) + 1/(2\gamma)||y - z||^2 }
#     = argmin_y { inf_x { f(x) : Lx = y } + 1/(2\gamma)||y - z||^2 }
#     = L * argmin_x { f(x) + 1/(2\gamma)||Lx - z||^2 }.
#
# When L is such that L'*L = mu*Id, then this just requires prox_{\gamma f}.
#
# In some other cases the prox can be "easily" computed, such as when f is
# quadratic or extended-quadratic.
#

export Epicompose

"""
**Epi-composition**

    Epicompose(L, f, [mu])

Returns the epi-composition of ``f`` with ``L``, also known as infimal
postcomposition or image function. Given a function f and a linear operator L,
their epi-composition is:
```math
g(y) = (Lf)(y) = inf_x { f(x) : Lx = y }
```
This is the dual operation to precomposition, see Rockafellar and Wets,
"Variational Analysis", Theorem 11.23.

If ``mu > 0`` is specified, then ``L`` is assumed to be such that ``L'*L == mu*I``,
and the proximal operator is computable for any convex ``f``. If ``mu`` is
not specified, then ``f`` must be of ``Quadratic`` type.
"""
abstract type Epicompose <: ProximableFunction end

Epicompose(L, f::T, mu) where {T <: ProximableFunction} = EpicomposeGramDiagonal(L, f, mu)
Epicompose(L, f::T) where {T <: Quadratic} = EpicomposeQuadratic(L, f)

# TODO add properties

### INCLUDE CONCRETE TYPES

include("epicomposeGramDiagonal.jl")
include("epicomposeQuadratic.jl")


text:
# Graph
mutable struct Graph{T} <: AbstractNexus{T}
    nn::Int
    na::Int
    fadj::Dict{T,Dict{T,Dict{Int,Nothing}}}
end

is_directed(::Graph{T}) where {T} = false
Graph{T}() where {T} = begin
    fadj = Dict{T,Dict{T,Dict{Int,Nothing}}}()
    return Graph(0, 0, fadj)
end

function empty!(g::Graph{T}) where {T}
    g.nn = 0
    g.na = 0
    empty!(g.fadj)
    return nothing
end

function add_arc!(G::Graph{T}, u::T, v::T, key::Union{Int,Nothing}=nothing) where {T}
    haskey(G.fadj, u) || begin
        G.fadj[u] = Dict()
        G.nn += 1
    end
    haskey(G.fadj, v) || begin
        G.fadj[v] = Dict()
        G.nn += 1
    end
    isnothing(key) && (key = new_arc_key(G, u, v))

    get!(G.fadj[u], v, Dict())
    get!(G.fadj[v], u, Dict())

    (haskey(G.fadj[u][v], key) || haskey(G.fadj[v][u], key)) || (G.na += 1)
    get!(G.fadj[u][v], key, nothing)
    get!(G.fadj[v][u], key, nothing)
    return key
end

# add_arc buat graph menambahkan u--v dan v--u ke fadj

function new_arc_key(g::Graph{T}, u::T, v::T) where {T}
    (haskey(g.fadj[u], v) || haskey(g.fadj[v], u)) ? begin
        keylist = union(keys(g.fadj[u][v]), keys(g.fadj[v][u]))
    end : return 1
    key = length(keylist)
    while key in keylist
        key += 1
    end # find lowest integer not used
    return key
end

function rem_arc!(G::Graph{T}, u::T, v::T, key::Union{Int,Nothing}=nothing) where {T}
    (haskey(G.fadj[u], v) || haskey(G.fadj[u], v)) ? begin
        d1 = G.fadj[u][v]
        d2 = G.fadj[v][u]
    end : return false # no v on u
    isnothing(key) ? begin
        pop!(d1)
        pop!(d2)
    end : (haskey(d1, key) || haskey(d2, key)) ? begin
        delete!(d1, key)
        delete!(d2, key)
    end : return false # rem key
    (length(d1) == 0 || length(d2) == 0) && begin
        delete!(G.fadj[u], v)
        delete!(G.badj[v], u)
    end
    G.na -= 1
    return true
end

function add_node!(g::Graph{T}, u) where {T}
    !haskey(g.fadj, u) ? begin # doesn't have node, create new dict
        g.fadj[u] = Dict()
        g.nn += 1
        return u
    end : return nothing
end

function rem_node!(g::Graph{T}, u::T) where {T}
    haskey(g.fadj, u) || return false

    for v in keys(g.fadj[u])
        g.na -= length(g.fadj[u][v])
        delete!(g.fadj[v], u)
        delete!(g.fadj[u], v)
    end
    delete!(g.fadj, u)
    g.nn -= 1
    return true
end


text:
@testset "COMMON OPTS" begin
	# -------------------- Test common_options ----------------------------------------
	@test GMT.get_GMTver()[1] >= v"6.0"
	GMT.__init__(true);
	GMT.dict2nt(Dict(:a =>1, :b => 2))
	@test GMT.parse_R(Dict(:xlim => (1,2), :ylim => (3,4), :zlim => (5,6)), "")[1] == " -R1/2/3/4/5/6"
	G1 = gmt("grdmath -R-2/2/-2/2 -I0.5 X Y MUL");
	@test GMT.build_opt_R(G1) == " -R-2/2/-2/2"
	@test GMT.build_opt_R(:d) == " -Rd"
	@test GMT.build_opt_R([]) == ""
	@test GMT.build_opt_R((bb=:global,)) == " -R-180/180/-90/90"
	@test GMT.build_opt_R((bb=:global360,)) == " -R0/360/-90/90"
	@test GMT.build_opt_R((bb=(1,2,3,4),)) == " -R1/2/3/4"
	@test GMT.build_opt_R((bb=(1,2,3,4), diag=1)) == " -R1/3/2/4+r"
	@test GMT.build_opt_R((bb_diag=(1,2,3,4),)) == " -R1/3/2/4+r"
	@test GMT.build_opt_R((continent=:s,)) == " -R=SA"
	@test GMT.build_opt_R((continent=:s,extend=4)) == " -R=SA+R4"
	@test GMT.build_opt_R((iso="PT,ES",extend=4)) == " -RPT,ES+R4"
	@test GMT.build_opt_R((iso="PT,ES",extend=[2,3])) == " -RPT,ES+R2/3"
	@test GMT.build_opt_R((bb=:d,unit=:k)) == " -Rd+uk"			# Idiot but ok
	@test_throws ErrorException("argument to the ISO key must be a string with country codes") GMT.build_opt_R((iso=:PT,))
	@test_throws ErrorException("No, no, no. Nothing useful in the region named tuple arguments") GMT.build_opt_R((zz=:x,))
	@test_throws ErrorException("Unknown continent name") GMT.build_opt_R((continent='a',extend=4))
	@test_throws ErrorException("Increments for limits must be a String, a Number, Array or Tuple") GMT.build_opt_R((iso="PT",extend='8'))
	@test GMT.opt_R2num(" -R1/2/3/4") == [1.0, 2, 3, 4]
	@test_throws ErrorException("The only valid case to provide a number to the 'proj' option is when that number is an EPSG code, but this (1500) is clearly an invalid EPSG")  GMT.build_opt_J(1500)
	@test GMT.build_opt_J(:X5)[1] == " -JX5"
	@test GMT.build_opt_J(2500)[1] == " -J2500"
	@test GMT.build_opt_J([])[1] == " -J"
	@test GMT.arg2str((1,2,3)) == "1/2/3"
	@test GMT.arg2str(("aa",2,3)) == "aa/2/3"
	@test GMT.arg2str(Dict(:shaded => "-4p/-6p/grey20@40"), [:shaded]) == "-4p/-6p/grey20@40"
	@test GMT.arg2str(Dict(:shaded => "aa bb"), [:shaded]) == "\"aa bb\""
	@test_throws ErrorException("arg2str: argument 'arg' can only be a String, Symbol, Number, Array or a Tuple, but was DataType") GMT.arg2str(typeof(1))
	@test GMT.parse_b(Dict(:b => (ncols=2, swapp_bytes=true, little_endian=true, type=:char)), "")[2] == " -b2cw+l"
	@test GMT.parse_b(Dict(:b => (ncols=2, type=:double)), "")[2] == " -b2d"
	@test GMT.parse_b(Dict(:b => (ncols=2, type=:ai)), "")[2] == " -b2d"
	@test GMT.parse_bi(Dict(:bi => (ncols=2, type=:float)), "")[2] == " -bi2f"
	@test GMT.parse_bo(Dict(:bo => (ncols=2, type=:int16)), "")[2] == " -bo2h"
	@test GMT.parse_c(Dict(:c => (1,2)), "")[1] == " -c0,1"
	@test GMT.parse_c(Dict(:c => [1,2]), "")[1] == " -c0,1"
	@test GMT.parse_c(Dict(:c => "1,2"), "")[1] == " -c0,1"
	@test GMT.parse_c(Dict(:c => 1), "")[1] == " -c0"
	@test GMT.parse_c(Dict(:c => "1"), "")[1] == " -c0"
	@test GMT.parse_l(Dict(:l => "ai ai"), "")[2] == " -l\"ai ai\""
	@test GMT.parse_l(Dict(:l => (text="ai ai", vspace=3)), "")[2] == " -l\"ai ai\"+G3"
	@test GMT.parse_n(Dict(:n => (bicubic=true,antialiasing=true,bc=:g)), "")[2] == " -nc+bg"
	@test GMT.parse_I(Dict(:inc => (x=1.5, y=2.6, unit="meter")),"", [:I :inc], "I") == " -I1.5e/2.6e"
	@test GMT.parse_I(Dict(:inc => (x=1.5, y=2.6, unit="m")),"", [:I :inc], "I") == " -I1.5m/2.6m"
	@test GMT.parse_I(Dict(:inc => (x=1.5, y=2.6, unit="data")),"", [:I :inc], "I") == " -I1.5/2.6u"
	@test GMT.parse_I(Dict(:inc => (x=1.5, y=2.6, extend="data")),"", [:I :inc], "I") == " -I1.5+e/2.6+e"
	@test GMT.parse_I(Dict(:inc => (x=1.5, y=2.6, unit="nodes")),"", [:I :inc], "I") == " -I1.5+n/2.6+n"
	@test GMT.parse_I(Dict(:inc => (2,4)),"", [:I :inc], "I") == " -I2/4"
	@test GMT.parse_I(Dict(:inc => [2 4]),"", [:I :inc], "I") == " -I2/4"
	@test GMT.parse_I(Dict(:inc => "2"),"", [:I :inc], "I") == " -I2"
	@test GMT.parse_I(Dict(:inc => "2"),"", [:I :inc], "") == "2"
	@test GMT.parse_JZ(Dict(:JZ => "5c"), "")[1] == " -JZ5c"
	@test GMT.parse_JZ(Dict(:Jz => "5c"), "")[1] == " -Jz5c"
	@test GMT.parse_JZ(Dict(:aspect3 => 1), " -JX10")[1] == " -JX10 -JZ10"
	@test GMT.parse_J(Dict(:J => "X5"), "", "", false)[1] == " -JX5"
	@test GMT.parse_J(Dict(:a => ""), "", "", true, true)[1] == " -J"
	@test GMT.parse_J(Dict(:J => "X", :figsize => 10), "")[1] == " -JX10"
	@test GMT.parse_J(Dict(:J => "X", :scale => "1:10"), "")[1] == " -Jx1:10"
	@test GMT.parse_J(Dict(:proj => "Ks0/15"), "")[1] == " -JKs0/15"
	@test GMT.parse_J(Dict(:scale=>"1:10"), "")[1] == " -Jx1:10"
	@test GMT.parse_J(Dict(:s=>"1:10"), "", " -JU")[1] == " -JU"
	@test GMT.parse_J(Dict(:J => "Merc", :figsize => 10), "", "", true, true)[1] == " -JM10"
	@test GMT.parse_J(Dict(:J => "+proj=merc"), "")[1] == " -J+proj=merc+width=" * split(GMT.def_fig_size, '/')[1]
	@test GMT.parse_J(Dict(:J => (name=:albers, parallels=[45 65])), "", "", false)[1] == " -JB0/0/45/65"
	@test GMT.parse_J(Dict(:J => (name=:albers, center=[10 20], parallels=[45 65])), "", "", false)[1] == " -JB10/20/45/65"
	@test GMT.parse_J(Dict(:J => "winkel"), "", "", false)[1] == " -JR"
	@test GMT.parse_J(Dict(:J => "M0/0"), "", "", false)[1] == " -JM0/0"
	@test GMT.parse_J(Dict(:J => (name=:merc,center=10)), "","", false)[1] == " -JM10"
	@test GMT.parse_J(Dict(:J => (name=:merc,parallels=10)), "","", false)[1] == " -JM0/0/10"
	@test GMT.parse_J(Dict(:J => (name=:Cyl_,center=(0,45))), "", "", false)[1] == " -JCyl_stere/0/45"
	@test_throws ErrorException("When projection arguments are in a NamedTuple the projection 'name' keyword is madatory.") GMT.parse_J(Dict(:J => (parallels=[45 65],)), "", "", false)
	@test_throws ErrorException("When projection is a named tuple you need to specify also 'center' and|or 'parallels'") GMT.parse_J(Dict(:J => (name=:merc,)), "", "", false)
	r = GMT.parse_params(Dict(:par => (MAP_FRAME_WIDTH=0.2, IO=:lixo, OI="xoli")), "");
	@test r == " --MAP_FRAME_WIDTH=0.2 --IO=lixo --OI=xoli"
	@test GMT.parse_params(Dict(:par => (:MAP_FRAME_WIDTH,0.2)), "") == " --MAP_FRAME_WIDTH=0.2"
	@test GMT.parse_params(Dict(:par => ("MAP_FRAME_WIDTH",0.2)), "") == " --MAP_FRAME_WIDTH=0.2"
	@test GMT.opt_pen(Dict(:lw => 5, :lc => :red),'W', nothing) == " -W5,red"
	@test GMT.opt_pen(Dict(:lw => 5),'W', nothing) == " -W5"
	@test GMT.opt_pen(Dict(:a => (10,:red)),'W', [:a]) == " -W10,red"
	@test_throws ErrorException("Nonsense in W option") GMT.opt_pen(Dict(:a => [1 2]),'W', [:a])
	@test GMT.get_color(((1,2,3),)) == "1/2/3"
	@test GMT.get_color(((1,2,3),100)) == "1/2/3,100"
	@test GMT.get_color(((0.1,0.2,0.35),)) == "26/51/89"
	@test GMT.get_color([1 2 3]) == "1/2/3"
	@test GMT.get_color([0.4 0.5 0.8; 0.1 0.2 0.75]) == "102/128/204,26/51/191"
	@test GMT.get_color([1 2 3; 3 4 5; 6 7 8]) == "1/2/3,3/4/5,6/7/8"
	@test GMT.get_color(:red) == "red"
	@test GMT.get_color((:red,:blue)) == "red,blue"
	@test GMT.get_color((200,300)) == "200,300"
	@test_throws ErrorException("GOT_COLOR, got an unsupported data type") GMT.get_color([1,2])
	@test_throws ErrorException("Color tuples must have only one or three elements") GMT.get_color(((0.2,0.3),))
	@test GMT.parse_unit_unit("data") == "u"
	@test GMT.parse_units((2,:p)) == "2p"
	@test GMT.add_opt((a=(1,0.5),b=2), (a="+a",b="-b")) == "+a1/0.5-b2"
	@test GMT.add_opt((symb=:circle, size=7, unit=:point), (symb="1", size="", unit="1")) == "c7p"
	@test GMT.add_opt(Dict(:L => "pen"), "", "L", [:L], (pen="_+p",)) == " -L+p"
	@test GMT.add_opt_1char("", Dict(:N=>"abc"), [[:N :geod2aux]]) == " -Na"
	@test GMT.add_opt_1char("", Dict(:N => ("abc", "sw", "x"), :Q=>"datum"), [[:N :geod2aux], [:Q :list]]) == " -Nasx -Qd"
	r = GMT.add_opt_fill("", Dict(:G=>(inv_pattern=12,fg="white",bg=[1,2,3], dpi=10) ), [:G :fill], 'G');
	@test r == " -GP12+b1/2/3+fwhite+r10"
	@test GMT.add_opt_fill("", Dict(:G=>:red), [:G :fill], 'G') == " -Gred"
	@test_throws ErrorException("For 'fill' option as a NamedTuple, you MUST provide a 'patern' member") GMT.add_opt_fill("", Dict(:G=>(inv_pat=12,fg="white")), [:G], 'G')
	d = Dict(:offset=>5, :bezier=>true, :cline=>"", :ctext=>true, :pen=>("10p",:red,:dashed));
	@test GMT.add_opt_pen(d, [:W :pen], "W") == " -W10p,red,dashed+cl+cf+s+o5"
	d = Dict(:W=>(offset=5, bezier=true, cline="", ctext=true, pen=("10p",:red,:dashed), arrow=(lenght=0.1,)));
	@test GMT.add_opt_pen(d, [:W :pen], "W") == " -W10p,red,dashed+cl+cf+s+o5+v"
	GMT.add_opt_cpt(Dict(:a=>1), "", [:b], 'A', 0, nothing, nothing, false, true, "-T0/10/1");
	GMT.add_opt_cpt(Dict(:a=>:red), "", [:a], 'A', 0, nothing, nothing, false, true, "-T0/10/1");

	r = vector_attrib(len=2.2,stop=[],norm="0.25i",shape=:arrow,half_arrow=:right,
	                  justify=:end,fill=:none,trim=0.1,endpoint=true,uv=6.6);
	@test r == "2.2+e+je+r+g-+n0.25i+h1+t0.1+s+z6.6"

	r = decorated(dist=("0.4i",0.25), symbol=:arcuate, pen=2, offset="10i", right=1);
	@test r == " -Sf0.4i/0.25+r+S+o10i+p2"
	r = decorated(dist=("0.8i","0.1i"), symbol=:star, symbsize=1, pen=(0.5,:green), fill=:blue, n_data=20, nudge=1, debug=1, dec2=1);
	@test r == " -S~d0.8i/0.1i:+sa1+d+gblue+n1+w20+p0.5,green"
	r = decorated(n_symbols=5, symbol=:star, symbsize=1, pen=(0.5,:green), fill=:blue, quoted=1);
	@test r == " -Sqn5:+p0.5,green"
	GMT.decorated((symbol="aiai",))		# Trigger a warning

	r = decorated(dist=("0.4i",0.25), angle=7, clearance=(2,3), debug=1, delay=1, font=10, color=:red, justify=:TC, const_label=:Ai, pen=(0.5,:red), fill=:blue, nudge=(3,4), rounded=1, unit=:TT, min_rad=0.5, curved=1, n_data=20, prefix="Pre", suffices="a,b", label=(:map_dist,"d"), quoted=1)
	@test r == " -Sqd0.4i/0.25:+a7+d+c2/3+e+f10+gred+jTC+lAi+n3/4+o+r0.5+uTT+v+w20+=Pre+xa,b+LDd+p0.5,red"

	@test GMT.helper_decorated(Dict(:line => [1 2 3 4; 5 6 7 8]), true) == "l1/2/3/4,5/6/7/8"
	@test GMT.helper_decorated(Dict(:Line => "CT/CB"), true) == "LCT/CB"
	@test GMT.helper_decorated(Dict(:Line => (1,2,"CB")), true) == "L1/2/CB"
	@test GMT.helper_decorated(Dict(:N_labels => 5), true) == "N5"
	@test_throws ErrorException("DECORATED: 'line' option. When array, it must be an Mx4 one") GMT.helper_decorated(Dict(:line => [1 2 3]), true)
	@test_throws ErrorException("DECORATED: 'dist' (or 'distance') option. Unknown data type.") GMT.helper_decorated(Dict(:dist => (a=1,)))
	@test GMT.parse_quoted(Dict(:label => "aa"), "") == "+Laa"
	@test GMT.parse_quoted(Dict(:label => :header), "") == "+Lh"
	@test GMT.parse_quoted(Dict(:label => :input), "") == "+Lf"
	@test_throws ErrorException("Wrong content for the :label option. Must be only :header or :input") GMT.parse_quoted(Dict(:label => :x), "")
	@test_throws ErrorException("Wrong content for the :label option. Must be only :plot_dist or :map_dist") GMT.parse_quoted(Dict(:label => (:x,)), "")
	GMT.parse_quoted(Dict(:label => 1),"")		# Trigger a warning
	GMT.helper_arrows(Dict(:geovec => "bla"));
	GMT.helper_arrows(Dict(:vecmap => "bla"));

	@test GMT.font(("10p","Times", :red)) == "10p,Times,red"
	r = text(text_record([0 0], "TopLeft"), R="1/10/1/10", J=:X10, F=(region_justify=:MC,font=("10p","Times", :red)), Vd=dbg2);
	ind = findfirst("-F", r); @test GMT.strtok(r[ind[1]:end])[1] == "-F+cMC+f10p,Times,red"

	@test GMT.build_pen(Dict(:lw => 1, :lc => [1,2,3])) == "1,1/2/3"
	@test GMT.parse_pen((0.5, [1 2 3])) == "0.5,1/2/3"

	@test GMT.helper0_axes((:left_full, :bot_full, :right_ticks, :top_bare, :up_bare)) == "WSetu"
	d=Dict(:xaxis => (axes=:WSen,title=:aiai, label=:ai, annot=:auto, ticks=[], grid=10, annot_unit=:ISOweek,seclabel=:BlaBla), :xaxis2=>(annot=5,ticks=1), :yaxis=>(custom="lixo.txt",), :yaxis2=>(annot=2,));
	@test GMT.parse_B(d, "")[1] == " -Bsya2 -Bsxa5f1 -Bpyclixo.txt -BpxaUfg10 -Bpx+lai+sBlaBla -BWSen+taiai"
	@test GMT.parse_B(Dict(:B=>:same), "")[1] == " -B"
	@test GMT.parse_B(Dict(:title => :bla), "")[1] == " -Baf -BWSen+tbla"
	@test GMT.parse_B(Dict(:frame => :auto, :title => :bla), "")[1] == " -Baf -BWSen+tbla"
	@test GMT.parse_B(Dict(:B=>:WS), "")[1] == " -Baf -BWS"
	@test GMT.parse_B(Dict(:title => "bla"), "", " -Baf")[1] == " -Baf -B+tbla"
	@test GMT.parse_B(Dict(:frame => (annot=10, title="Ai Ai"), :grid => (pen=2, x=10, y=20)), "", " -Baf -BWSen")[1] == " -Bpa10 -Byg20 -Bxg10 -BWSen+t\"Ai Ai\""
	@test GMT.parse_B(Dict(:frame => (axes=(:left_full, :bottom_full, :right_full, :top_full), annot=10)), "")[1] == " -Bpa10 -BWSEN"
	@test GMT.parse_B(Dict(:xaxis => (axes=:full, annot=10)), "")[1] == " -Bpxa10 -BWSEN"
	@test GMT.parse_B(Dict(:frame => (fill=220,)), "", " -Baf -Bg -BWSne")[1] == " -Baf -Bg -BWSne+g220"
	@test GMT.parse_B(Dict(:frame => :full), "")[1] == " -Baf -BWSEN"
	GMT.helper2_axes("lolo");
	@test_throws ErrorException("Custom annotations NamedTuple must contain the member 'pos'") GMT.helper3_axes((a=0,),"","")

	@test GMT.consolidate_Baxes(" -Baf -BWSen -BpxaUfg10 -BWSen+taiai -Bpx+lai+sBlaBla -Bpyclixo.txt -Bsxa5f1") ==
		" -BWSen -BpxaUfg10 -BWSen+taiai -Bpx+lai+sBlaBla -Bpyclixo.txt -Bsxa5f1"
	@test GMT.consolidate_Baxes(" -Baf -Bza -BWSrt -Bpxapi2f0.5 -Bpx+lpis -Bpya1 -Bpy+a60") ==
		" -Byf -Bza -BWSrt -Bpxapi2f0.5+lpis -Bpya1+a60"
	@test GMT.helper_consolidate_B("","a",false,false) == ""
	@test GMT.consolidate_Bframe(" -Bpx+lx -B+gwhite -Baf -BWSen -By+lx") == " -Bpx+lx -Baf -By+lx -BWSen+gwhite"
	@test GMT.consolidate_Bframe(" -Bpx+lx -Bpy+lx -BWSrt+gwhite") == " -Bpx+lx -Bpy+lx -BWSrt+gwhite"
	@test GMT.consolidate_Bframe(" -Bpx+lx -BWSrt+gwhite -Bpy+lx") == " -Bpx+lx -Bpy+lx -BWSrt+gwhite"

	d=Dict(:L => (pen=(lw=10,lc=:red),) );
	@test GMT.add_opt(d, "", "", [:L], (pen=("+p",GMT.add_opt_pen),) ) == "+p10,red"
	r = psxy([0.0, 1],[0, 1.1], L=(pen=(10,:red),bot=true), Vd=dbg2);
	@test startswith(r,"psxy  -JX" * GMT.def_fig_size * " -Baf -BWSen -R-0.04/1.04/-0.04/1.12 -L+p10,red+yb")
	r = psxy([0.0, 1],[0, 1.1], L=(pen=(lw=10,cline=true),bot=true), Vd=dbg2);
	@test startswith(r,"psxy  -JX" * GMT.def_fig_size * " -Baf -BWSen -R-0.04/1.04/-0.04/1.12 -L+p10+cl+yb")
	@test startswith(psxy([0.0, 1],[0, 1.1], figsize=(10,12), aspect=:equal, Vd=dbg2), "psxy  -JX10/12")
	@test startswith(psxy([0.0, 1],[0, 1.1], figsize=10, aspect=:equal, Vd=dbg2), "psxy  -JX10/0")
	@test startswith(psxy([0.0, 1],[0, 1.1], aspect=:equal, Vd=dbg2), "psxy  -JX" * split(GMT.def_fig_size, '/')[1] * "/0")
	psxy!([0 0; 1 1.1], Vd=dbg2);
	psxy!("", [0 0; 1 1.1], Vd=dbg2);
	GMT.get_marker_name(Dict(:y => "y"), nothing, [:y], false)
	@test GMT.line_decorated_with_symbol(Dict()) == " -S~d0.88:+sc0.11+gwhite+p0.75,black"
	@test_throws ErrorException("Argument of the *bar* keyword can be only a string or a NamedTuple.") GMT.parse_bar_cmd(Dict(:a => 0), :a, "", "")

	@test_throws ErrorException("Custom annotations NamedTuple must contain the member 'pos'") GMT.helper3_axes((post=1:5,), "p", "x")
	GMT.helper3_axes(1,"","")		# Trigger a warning

	dt = collect(Dates.DateTime(Dates.now()):Dates.Month(6):Dates.DateTime(Dates.now() + Dates.Year(10)));
	GMT.read_data(Dict{Symbol, Any}(), "", "", dt)
	GMT.read_data(Dict{Symbol, Any}(), "", "", [dt rand(length(dt))])

	GMT.round_wesn([1.333 17.4678 6.66777 33.333], true);
	GMT.round_wesn([1 1 2 2]);
	GMT.round_wesn([1. 350 1. 180], true)
	GMT.round_wesn([0. 1.1 0. 0.1], true)

	GMT.GMTdataset([0.0 0]);
	GMT.GMTdataset([0.0 0], Array{String,1}());

	img16 = rand(UInt16, 16, 16, 3);
	I = GMT.mat2img(img16);
	img16 = rand(UInt16, 4, 4, 3);
	I = GMT.GMTimage("", "", 0, [1.,4,1,4,0,255], [1., 1], 1, NaN, "", String[], String[], collect(1.:4),collect(1.:4),zeros(3),img16, vec(zeros(Clong,1,3)), 0, Array{UInt8,2}(undef,1,1), "TCBa", 0)
	GMT.mat2img(I);
	GMT.mat2img(img16, histo_bounds=8440);
	GMT.mat2img(img16, histo_bounds=[8440 13540]);
	GMT.mat2img(img16, histo_bounds=[8440 13540 800 20000 1000 30000]);
	GMT.mat2img(rand(UInt16,32,32,3),stretch=:auto);
	I = GMT.mat2img(img16, I)
	slicecube(I, 2)
	@test rescale(1:5,-1,1) == -1.0:0.5:1
	@test rescale(1:4, type=UInt8) == [0, 85, 170, 255]
	magic(4)
	magic(6)

	D = mat2ds([0 0; 1 1],["a", "b"]);	D.header = "a";
	GMT.make_zvals_vec(D, ["a", "b"], [1,2]);
	GMT.edit_segment_headers!([D], [1], "0")
	GMT.edit_segment_headers!(D, [1], "0")

	@test_throws ErrorException("Bad 'stretch' argument. It must be a 1, 2 or 6 elements array and not 3") GMT.mat2img(img16, histo_bounds=[8440 13540 0]);
	@test_throws ErrorException("Memory layout option must have 3 characters and not 1") GMT.parse_mem_layouts("-%1")
	@test_throws ErrorException("Memory layout option must have at least 2 chars and not 1") GMT.parse_mem_layouts("-&1")
	@test GMT.parse_mem_layouts("-&BR")[3] == "BR"
	@test_throws ErrorException("parse_arg_and_pen: Nonsense first argument") GMT.parse_arg_and_pen(([:a],0))
	@test_throws ErrorException("GMT: No module by that name -- bla -- was found.") gmt("bla")
	#@test_throws ErrorException("grd_init: input (Int64) is not a GRID container type") GMT.grid_init(C_NULL,0,0)
	@test_throws ErrorException("image_init: input is not a IMAGE container type") GMT.image_init(C_NULL,0)
	@test_throws ErrorException("Bad family type") GMT.GMT_Alloc_Segment(C_NULL, -1, 0, 0, "", C_NULL)
	#@test_throws ErrorException("Unknown family type") GMT.GMT_Create_Data(C_NULL, -99, 0, 0)
	@test_throws ErrorException("Expected a PS structure for input") GMT.ps_init(C_NULL, 0, 0)
	@test_throws ErrorException("size of x,y vectors incompatible with 2D array size") GMT.grdimg_hdr_xy(rand(3,3), 0, 0, [1 2], [1])
	GMT.strncmp("abcd", "ab", 2)
	GMT.parse_proj((name="blabla",center=(0,0)))

	@test GMT.parse_i(Dict(:i=>(0,1,2,2)), "")[1] == " -i0,1,2,2"
	@test GMT.parse_j(Dict(:spherical_dist => "f"), "")[1] == " -jf"
	@test GMT.parse_t(Dict(:t=>0.2), "")[1] == " -t20.0"
	@test GMT.parse_t(Dict(:t=>20), "")[1]  == " -t20"
	@test GMT.parse_contour_AGTW(Dict(:A => [1]), "")[1] == " -A1,"
	GMT.helper2_axes("");
	@test GMT.axis(ylabel="bla")[1] == " -Bpy+lbla";
	@test GMT.axis(Yhlabel="bla")[1] == " -Bpy+Lbla";
	@test GMT.axis(scale="exp")[1] == " -Bpp";
	@test GMT.axis(phase_add=10)[1] == " -Bp+10";
	@test GMT.axis(phase_sub=10)[1] == " -Bp-10";

	@test GMT.parse_grid(Dict(), (z=25,)) == " -Bzg25"
	@test GMT.parse_grid(Dict(), ()) == " -Bg"
	@test GMT.parse_grid(Dict(), (xyz=1,)) == " -Bg -Bzg"
	@test GMT.parse_grid(Dict(), "x") == " -Bxg"
	@test GMT.parse_grid(Dict(), "x25") == " -Bxg25"
	@test GMT.parse_grid(Dict(), "y25") == " -Byg25"
	@test GMT.parse_grid(Dict(), "xyz") == " -Bg -Bzg"
	@test GMT.parse_grid(Dict(), :on) == " -Bg"
	r = coast(region=:global, proj=:Winkel, par=(MAP=1,), grid=(pen=:red,), Vd=2)
	@test occursin("--MAP=1 --MAP_GRID_PEN_PRIMARY=red", r)

	r,o = GMT.prepare2geotif(Dict(:geotif => :trans), ["pscoast  -Rd -JX12cd/0 -Baf -BWSen -W0.5p -Da"], "", false);
	@test startswith(r[1],"pscoast  -Rd -JX30cd/0 -W0.5p -Da  -B0 --MAP_FRAME_TYPE=inside --MAP_FRAME_PEN=0.1,254")
	@test o == " -TG -W+g"
	r,o = GMT.prepare2geotif(Dict(:kml => :trans), ["pscoast  -Rd -JX12cd/0 -Baf -BWSen -W0.5p -Da"], "", false);
	@test startswith(r[1],"pscoast  -Rd -JX30cd/0 -W0.5p -Da  -B0 --MAP_FRAME_TYPE=inside --MAP_FRAME_PEN=0.1,254")
	@test o == " -TG -W+k"
	o = GMT.prepare2geotif(Dict(:kml => "+tLoLo"), ["pscoast  -Rd -JX12cd/0 -Baf -W0.5p -Da"], "", false)[2];
	@test o == " -TG -W+k+tLoLo"
	o = GMT.prepare2geotif(Dict(:kml => (title=:Lolo, layer=:bla, fade=(1,2), URL="http")), ["pscoast  -Rd -JX12 -Baf -W0.5p -Da"], "", false)[2];
	@test o == " -TG -W+k+tLolo+nbla+f1/2+uhttp"
	coast(region=:global, kml=:trans, proj=:merc,Vd=dbg2)

	@test (GMT.check_axesswap(Dict(:axesswap => (y=1,)), "?") == "-?")
	@test (GMT.check_axesswap(Dict(:axesswap => ("x", :y)), "?/?") == "-?/-?")
	@test (GMT.check_axesswap(Dict(:axesswap => (xy=1,)), "?/?") == "-?/-?")
	@test (GMT.check_axesswap(Dict(:axesswap => ("xy")), "?/?") == "-?/-?")
	@test (GMT.check_axesswap(Dict(:axesswap => ("y")), "?/?") == "?/-?")

	@test (GMT.parse_opt_range(Dict(:T => (1,2,0.1,:num)), "", "")[1] == "1/2/0.1+n")
	@test (GMT.parse_opt_range(Dict(:T => (1,2,0.1,:num)), "", "T")[1] == " -T1/2/0.1+n")
	@test (GMT.parse_opt_range(Dict(:T => (1,2,0.1,:log1)), "")[1] == "1/2/0.1+l")
	@test (GMT.parse_opt_range(Dict(:T => [1]), "")[1] == "1,")
	GMT.parse_opt_range(Dict(:T => (1,2,0.1,:mum,:log2)), "")[1]	# Prints a warning

	GMT.round_datetime([DateTime(2013,1,1), DateTime(2013,1,1,0,0,1)]);

	GMT.GMT_PEN();
	GMT.GMT_PEN(0.0, 0.0, (0.0, 0.0, 0.0, 0.0), map(UInt8, (repeat('\0', 128)...,)), 0, 0, (pointer([0]), pointer([0])));

	GMT.guess_proj([0.0 0.0], [0.0 0.0])
	GMT.guess_proj([0., 20.], [0.0, 20.])
	GMT.guess_proj([0., 20.], [35.0, 45.])
	GMT.guess_proj([0., 20.], [80.0, 90.])
	GMT.guess_proj([0., 20.], [-90.0, -80.])
	GMT.guess_proj([0., 20.], [-6.0, 90.])
	GMT.guess_proj([0., 20.], [-30.0, 30.])
	GMT.guess_proj([0., 20.], [-40.0, 60.])
	GMT.guess_WESN(Dict(:p=>(350,2)), "")
	GMT.guess_WESN(Dict(:p=>"350/2"), "")
	GMT.parse_q(Dict(:p=>(350,2)), "")

	@test_throws ErrorException("Only integer or floating point types allowed in input. Not this: Char") GMT.dataset_init(GMT.API, ' ', [0])
	@test_throws ErrorException("Wrong type (Int64) for the 'text' argin") GMT.text_record(rand(2,2), 0)

	GMT.show_non_consumed(Dict(:lala => 0), "prog");
	GMT.dbg_print_cmd(Dict(:lala => 0, :Vd=>2), "prog");

	GMT.justify("aiai")		# A warning

	GMT.grid2pix(GMT.peaks())
	GMT.grid2pix([0., 1, 0, 1, 0, 1, 1, 1, 1], pix=false)

	gmthelp([:n :sphinterpolate])
	gmthelp(:wW)

	# Test here is to the showfig fun
	grdimage([1 2;3 4])
	showfig(savefig="lixo.png",show=false)

	y = [NaN 2 3 4;5 6 NaN 8;9 10 11 12];
	@test nanmean(y,1) == [7.0  6.0  7.0  8.0]
	nanstd(y,1)
	isnodata(rand(2,2))

	@test doy2date(283, 2021) == Date("2021-10-10")
	@test date2doy("2021-10-10") == 283
	@test GMT.yeardecimal("2000") == 2000.0

	# STACKGRIDS
	GMT.stackgrids(["a","b"], [1,2], mirone=true)
	rm("automatic_list.txt")
	gmtwrite("lixo1.grd", mat2grid(rand(Float32, 16, 16), proj4=GMT.prj4WGS84))
	gmtwrite("lixo2.grd", mat2grid(rand(Float32, 16, 16), proj4=GMT.prj4WGS84))
	GMT.stackgrids(["lixo1.grd", "lixo2.grd"], [now(), now()+Dates.Day(1)], save="lixo_cube.nc", z_unit="date")
	rm("lixo1.grd")
	rm("lixo2.grd")
	#rm("lixo_cube.nc")		# can't be deleted because Julia still holds its file handle.

	gmtbegin()
	GMT.gmt_restart(false)
	resetGMT()
end


text:
# Infer population model parameters from control group

using Pumas, MCMCChains, StatsPlots, CSV, Random
plotly()

# Import data
path = joinpath(
    @__DIR__,
    "../../..", "python/lung_cancer/data",
    "lxf_control_growth.csv")
data = DataFrame(CSV.File(path))

# Filter Body weight information
data = filter(x -> !(occursin.("Body weight", x.Biomarker)), data)
pumas_data = read_pumas(data,
    id           = :ID,
    time         = :Time,
    observations = [:Measurement],
    event_data = false)

# Visualise tumour volume
@df data scatter(
    :Time,
    :Measurement,
    group = :ID)
yaxis!("Tumour volume in cm^3")
xaxis!("Time in day")

# Build model
tumour_growth_model = @model begin
    @param   begin
        mean_initial_volume ~ Constrained(
            Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        mean_critical_volume ~ Constrained(
            Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        mean_growth_rate ~ Constrained(
            Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        std ~ Constrained(
            MvNormal([1 , 1, 1]),
            lower=[0, 0, 0],
            upper=[Inf, Inf, Inf],
            init=[0.09,0.09, 0.09])
        σ_base ~ Constrained(
            Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        σ_rel ~ Constrained(
            Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        # mean_initial_volume ∈ RealDomain(lower=0, init = 0.2)
        # mean_critical_volume ∈ RealDomain(lower=0, init = 1)
        # mean_growth_rate ∈ RealDomain(lower = 0, init= 0.1)
        # std ∈ PDiagDomain(init=[0.09,0.09, 0.09])
        # σ_base ∈ RealDomain(lower=0, init=0.04)
        # σ_rel ∈ RealDomain(lower=0, init=0.04)
    end

    @random begin
        η ~ MvNormal(std)
    end

    @pre begin
        initial_volume = mean_initial_volume * exp(η[1])
        critical_volume = mean_critical_volume * exp(η[2])
        growth_rate = mean_growth_rate * exp(η[3])
    end

    @init begin
        volume = initial_volume
    end

    @dynamics begin
       volume' = (growth_rate * critical_volume * volume) / (volume + critical_volume)
    end

  @derived begin
      σ_tot := σ_base .+ volume * σ_rel
      Measurement ~ @. Normal(volume, σ_tot)
    end
end

# Infer parameters
param = init_param(tumour_growth_model)
result = fit(tumour_growth_model, pumas_data, param, Pumas.BayesMCMC();
  nsamples=2000, nadapts=1000)

# Visualise population parameters
chains = Pumas.Chains(result)
plot(chains[1000:end])

max_stars_repo_name: DavAug/erlotinib

max_issues_repo_name: DavAug/erlotinib

max_forks_repo_name: DavAug/erlotinib


text:
"""
    card(deck, idx)

Return the card at position `idx` from the given `deck`.
"""
card(deck, idx) = deck[idx]

"""
    replace_card!(deck, idx => new_card)

Exchange the card at position `idx` with `new_card` and return the adjusted `deck`.
"""
function replace_card!(deck, r)
    deck[r.first] = r.second
    deck
end

"""
    insert_card_at_top!(deck, new_card)

Insert `new_card` at the top of `deck` and return the modified `deck`.
"""
insert_card_at_top!(deck, new_card) = push!(deck, new_card)

"""
    remove_card!(deck, idx)

Remove the card at position `idx` from `deck` and return the modified `deck`.
"""
remove_card!(deck, idx) = deleteat!(deck, idx)

"""
    remove_card_from_top!(deck)

Remove the card at the top of `deck` and return the modified `deck`.
"""
function remove_card_from_top!(deck)
    pop!(deck)
    deck
end

"""
    insert_card_at_bottom!(deck, new_card)

Insert `new_card` at the bottom of `deck` and return the modified `deck`.
"""
insert_card_at_bottom!(deck, new_card) = pushfirst!(deck, new_card)

"""
    remove_card_from_bottom!(deck)

Remove the card at the bottom of `deck` and return the modified `deck`.
"""
function remove_card_from_bottom!(deck)
    popfirst!(deck)
    deck
end

"""
    check_deck_size(deck, deck_size)

Check if the size of `deck` equals the given `deck_size`.
"""
check_deck_size(deck, deck_size) = length(deck) == deck_size


text:
"""
    QXZoo.jl - implementation of circuit generator for
    different quantum algorithms.
"""
module QXZoo

using Reexport

include("gates/gate_ops.jl")
include("gates/default_gates.jl")
include("gates/gate_map.jl")
include("gates/composite_gates.jl")

include("circuits/circuit.jl")
include("algorithms/Algorithms.jl")

@reexport using QXZoo.GateMap
@reexport using QXZoo.GateOps
@reexport using QXZoo.DefaultGates
@reexport using QXZoo.CompositeGates
@reexport using QXZoo.Circuit

# Initialise default gates into cache
GateMap.init_cache()

end


text:
"""

    air()

Produce an array of sample climate data for testing. This sample time
series is taken from the [NCEP/NCAR Reanalysis 1 dataset](https://doi.org/10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2),
representing the 925 hPa air temperature standardized anomalies of all
DJF winters from 1948-2018 at 50.0°N, 235.0°E.

See also:

>Kalnay, E., and Coauthors, 1996: The NCEP/NCAR 40-year reanalysis project. *Bulletin of the American Meteorological Society*, **77**, 437–471, [https://doi.org/10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2](https://doi.org/10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2).

# Examples
```jldoctest
julia> air()
6408-element Array{Float64,1}:
  0.67547214
  0.8308915
  0.19012427
 -0.20711477
 -0.10181273
  0.7421673
  0.22026125
  0.050515447
  0.043254483
  0.044615462
  ⋮
  0.7591071
  0.54696834
  0.3421846
  0.044188347
  0.2433472
  0.4460243
  1.4941149
  0.49388295
 -0.052520715

```

"""
function air()
    return vec(readdlm(joinpath(dirname(pathof(SGSDist)),
                                "../test/data/R1_air_stdanom_t1948-2018-djf_x235p0_y50p0_z925.csv"),' ','\n'))
end


text:
# example: julia opf.jl 66200 /path/to/JuMPSupplement/acpower
# JuMPSupplement is at https://github.com/mlubin/JuMPSupplement

#include("compare_jump.jl")

using JuMP
using Compat
using Ipopt
using ReverseDiffSparse2

type Bus
    bustype::Int
    name::ASCIIString
    voltage0::Float64
    angle0::Float64
    p_gen::Float64
    q_gen::Float64
    q_min::Float64
    q_max::Float64
    p_load::Float64
    q_load::Float64
    g_shunt::Float64
    b_shunt0::Float64
    b_shunt_min::Float64
    b_shunt_max::Float64
    b_dispatch::Float64
    area::Float64
end

type Branch
    from::Int
    to::Int
    branchtype::Int
    r::Float64
    x::Float64
    c::Float64
    tap0::Float64
    tap_min0::Float64
    tap_max0::Float64
    def0::Float64
    def_min::Float64
    def_max::Float64
    g::Float64
    b::Float64
end

function opf(bus, branch,
    bus_voltage_min, bus_voltage_max,
    branch_tap_min, branch_tap_max,
    p_gen_upper, p_gen_lower)

    mod = Model()

    nbus = length(bus)
    nbranch = length(branch)
    @defVar(mod, bus_voltage_min[bus[i].bustype] <= bus_voltage[i=1:nbus] <= bus_voltage_max[bus[i].bustype], start=1)
    @defVar(mod, bus[i].b_shunt_min <= bus_b_shunt[i=1:nbus] <= bus[i].b_shunt_max, start=bus[i].b_shunt0)
    @defVar(mod, bus_angle[1:nbus], start=0)

    @defVar(mod, branch_tap_min <= branch_tap[1:nbranch] <= branch_tap_max, start=1)
    @defVar(mod, branch[i].def_min <= branch_def[i=1:nbranch] <= branch[i].def_max, start = branch[i].def0)

    in_lines = [Int[] for i in 1:nbus] # indices of incoming branches
    out_lines = [Int[] for i in 1:nbus] # indices of outgoing branches
    @show nbranch
    @show nbus
    for i in 1:nbranch
        b = branch[i]
        push!(out_lines[b.from],i)
        push!(in_lines[b.to],i)
        @assert 1 <= b.to <= nbus
    end

    @defNLExpr(mod, Gself[k=1:nbus],
        bus[k].g_shunt +
        sum{branch[i].g*branch_tap[i]^2, i=out_lines[k]} +
        sum{ branch[i].g, i=in_lines[k]})

    @defNLExpr(mod, Gout[i=1:nbranch],
        (-branch[i].g*cos(branch_def[i])+branch[i].b*sin(branch_def[i]))*branch_tap[i])
    @defNLExpr(mod, Gin[i=1:nbranch],
        (-branch[i].g*cos(branch_def[i])-branch[i].b*sin(branch_def[i]))*branch_tap[i])

    @defNLExpr(mod, Bself[k=1:nbus],
        bus_b_shunt[k] +
        sum{branch[i].b*branch_tap[i]^2 + branch[i].c/2, i=out_lines[k]} +
        sum{ branch[i].b + branch[i].c/2, i=in_lines[k]})

    @defNLExpr(mod, Bin[i=1:nbranch],
        (branch[i].g*sin(branch_def[i])-branch[i].b*cos(branch_def[i]))*branch_tap[i])
    @defNLExpr(mod, Bout[i=1:nbranch],
        (-branch[i].g*sin(branch_def[i])-branch[i].b*cos(branch_def[i]))*branch_tap[i])

    # Minimize active power

    @setNLObjective(mod, Min, sum{ (bus[k].p_load +
        sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
            (Gin[i] * cos(bus_angle[k] - bus_angle[branch[i].from]) +
             Bin[i] * sin(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } +
        sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
            (Gout[i] * cos(bus_angle[k] - bus_angle[branch[i].to]) +
             Bout[i] * sin(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] }
        + bus_voltage[k]^2*Gself[k] )^2,
        k in 1:nbus; bus[k].bustype == 2 || bus[k].bustype == 3})
    for k in 1:nbus
        bus[k].bustype == 0 || continue
        # p_load
        @addNLConstraint(mod, bus[k].p_gen - bus[k].p_load -
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * cos(bus_angle[k] - bus_angle[branch[i].from]) +
                 Bin[i] * sin(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } -
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * cos(bus_angle[k] - bus_angle[branch[i].to]) +
                 Bout[i] * sin(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] }
            - bus_voltage[k]^2*Gself[k] == 0)
        # q_load
        @addNLConstraint(mod, bus[k].q_gen - bus[k].q_load -
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * sin(bus_angle[k] - bus_angle[branch[i].from]) -
                 Bin[i] * cos(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } -
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * sin(bus_angle[k] - bus_angle[branch[i].to]) -
                 Bout[i] * cos(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] }
             + bus_voltage[k]^2*Bself[k] == 0)
    end
    for k in 1:nbus
        (bus[k].bustype == 2 || bus[k].bustype == 3) || continue
        # q_inj
        @addNLConstraint(mod, bus[k].q_min <= bus[k].q_load +
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * sin(bus_angle[k] - bus_angle[branch[i].from]) -
                 Bin[i] * cos(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } +
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * sin(bus_angle[k] - bus_angle[branch[i].to]) -
                 Bout[i] * cos(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] } +
            - bus_voltage[k]^2*Bself[k] <= bus[k].q_max)
        # p_inj
        @addNLConstraint(mod, 0 <= bus[k].p_load +
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * cos(bus_angle[k] - bus_angle[branch[i].from]) +
                 Bin[i] * sin(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } +
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * cos(bus_angle[k] - bus_angle[branch[i].to]) +
                 Bout[i] * sin(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] } +
            + bus_voltage[k]^2*Gself[k] <= p_gen_upper*bus[k].p_gen)
    end

    for i in 1:nbus
        #---- FREEZE THE REFERENCE BUS ANGLE TO ZERO.
        if bus[i].bustype == 3
            setLower(bus_angle[i], 0.0)
            setUpper(bus_angle[i], 0.0)
        end
        #---- FREEZE ANY DISPATCHABLE SHUNTS.
        if bus[i].b_dispatch == 0
            setLower(bus_b_shunt[i],bus[i].b_shunt0)
            setUpper(bus_b_shunt[i],bus[i].b_shunt0)
        end
    end

    for i in 1:nbranch
        #---- FREEZE ANY BRANCH TAPS.
        if branch[i].branchtype == 0 || branch[i].branchtype == 3
            setLower(branch_tap[i],1)
            setUpper(branch_tap[i],1)
        end
        #---- FREEZE CERTAIN PHASE SHIFTERS.
        if branch[i].branchtype != 4
            setLower(branch_def[i],getValue(branch_def[i]))
            setUpper(branch_def[i],getValue(branch_def[i]))
        end
    end

    return mod

end

function prepdata(bus, branch)

    buses = Array(Bus,size(bus,1))
    refbus = 0
    shunt = 0
    busmap = Dict()
    for i in 1:size(bus,1)
        @assert !haskey(busmap, bus[i,1])
        busmap[bus[i,1]] = i
        buses[i] = Bus(bus[i,2:end]...)
        # rescale
        buses[i].p_gen /= 100
        buses[i].q_gen /= 100
        buses[i].q_min /= 100
        buses[i].q_max /= 100
        buses[i].p_load /= 100
        buses[i].q_load /= 100
        if buses[i].bustype == 3
            refbus += 1
        end
        if buses[i].b_dispatch == 0
            shunt += 1
        end
    end
    @show refbus
    @show shunt

    branches = Array(Branch,size(branch,1))
    branchtap = 0
    phaseshift = 0
    for i in 1:size(branch,1)
        branches[i] = Branch(branch[i,2:end]...,0,0)
        @assert branch[i,1] == i
        branches[i].to = busmap[branches[i].to]
        branches[i].from = busmap[branches[i].from]
        branches[i].g = branches[i].r /(branches[i].r^2 + branches[i].x^2)
        branches[i].b = -branches[i].x /(branches[i].r^2 + branches[i].x^2)
        # rescale
        branches[i].def_min *= 3.14159/180
        branches[i].def_max *= 3.14159/180
        branches[i].def0 *= -3.14159/180
        if branches[i].branchtype == 0 || branches[i].branchtype == 3
            branchtap += 1
        end
        if branches[i].branchtype != 4
            phaseshift += 1
        end
    end
    @show branchtap
    @show phaseshift
    return branches, buses
end

numbuses = ARGS[1]
path = ARGS[2]
branches, buses = prepdata(readdlm("$path/IEEE$numbuses.bus"),readdlm("$path/IEEE$numbuses.branch"))

bus_voltage_min = @compat Dict(0 => 0.85, 1 => 0.85, 2 => 0.92, 3 => 0.99)
bus_voltage_max = @compat Dict(0 => 1.15, 1 => 1.15, 2 => 1.08, 3 => 1.01)
branch_tap_min = 0.85
branch_tap_max = 1.15

p_gen_upper = 1.10
p_gen_lower = 0.90
mod = opf(buses, branches,
    bus_voltage_min, bus_voltage_max,
    branch_tap_min, branch_tap_max,
    p_gen_upper, p_gen_lower)

setSolver(mod,IpoptSolver(max_iter=100))
#setSolver(mod,RDSSolver(IpoptSolver()))
println("solving")
solve(mod)
#compare_jump(mod)


text:
# ---
# title: 113. Path Sum II
# id: problem113
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/path-sum-ii/description/>
# hidden: true
# ---
#
# Given a binary tree and a sum, find all root-to-leaf paths where each path's
# sum equals the given sum.
#
# **Note:**  A leaf is a node with no children.
#
# **Example:**
#
# Given the below binary tree and `sum = 22`,
#
#
#
#           **5**
#          **/ \**
#         **4   8**
#        **/**   / **\**
#       **11**  13  **4**
#      /  **\**    **/** \
#     7    **2**  **5**   1
#
#
# Return:
#
#
#
#     [
#        [5,4,11,2],
#        [5,8,4,5]
#     ]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    resolvedfield(d,u)
    Construct appropriately conjugated resolved field for any field type
"""
function resolvedfield(d::AbstractDomain,u::Array{T,N}) where {T,N}
    (nx,ny) = size(d)
    û = zeros(T,2ny-1,2nx-1)
    Λ = N == 2 ? size(u)[2] - 1 : 0 # interpret cutoff for DNS,DSS,GSS
    for m1=0:Λ
        nmin = m1==0 ? 1 : -ny+1
        for n1 = nmin:ny-1
            û[n1+ny,m1+nx] = u[n1+ny,m1+1]
            û[-n1+ny,-m1+nx] = conj(u[n1+ny,m1+1])
        end
    end
    û
end

resolvedfield(d::AbstractDomain,u::Union{DSSField{T},GSSField{T}}) where {T<:AbstractFloat} = resolvedfield(d,u.x[1])

"""
    zonalfield(d,u)
    Construct appropriately conjugated zonal component of any field type
"""
function zonalfield(d::AbstractDomain,u::Array{T,1}) where {T}
    ny = d.ny
    û = zeros(T,2ny-1)
    for n = 1:ny-1
        û[n+ny] = u[n+ny]
        û[-n+ny] = conj(u[n+ny])
    end
    û
end

zonalfield(d::AbstractDomain,u::DNSField{T}) where {T<:AbstractFloat} = zonalfield(d,u[:,1])
zonalfield(d::AbstractDomain,u::DSSField{T}) where {T<:AbstractFloat} = zonalfield(d,u.x[1])
zonalfield(d::AbstractDomain,u::GSSField{T}) where {T<:AbstractFloat} = zonalfield(d,u.x[1][:,1])

"""
    inversefourier(d,u)
    Inverse Fourier functions of appropriately conjugated fields
"""

function inversefourier(d::AbstractDomain,u::FirstCumulant{T}) where {T<:AbstractFloat}
    s = (2d.ny-1)/2.0
    s*real(ifft(ifftshift(u)))
end

function inversefourier(d::AbstractDomain,u::Field{T}) where {T<:AbstractFloat}
    s = (2d.ny-1)*(2d.nx-1)/4.0
    s*real(ifft(ifftshift(u)))
end

"""
    vorticity(d,u)
    zonalvorticity(d,u)
    Compute vorticity and zonal vorticity based on input solution
"""
vorticity(d,u) = resolvedfield(d,u) |> x->inversefourier(d,x)
zonalvorticity(d,u) = zonalfield(d,u) |> x->inversefourier(d,x) # the beauty of julia!

"""
    Velocity
"""
function xvelocityfield(d::AbstractDomain,u::Array{Complex{T},N}) where {T<:AbstractFloat,N}
    (lx,ly),(nx,ny) = length(d),size(d)
    û = fill!(similar(u),0)
    Λ = N == 2 ? size(u)[2] - 1 : 0 # interpret cutoff for DNS,DSS,GSS
    for m1=0:Λ
        nmin = m1==0 ? 1 : -ny+1
        for n1 = nmin:ny-1
            kx,ky = 2π*m1/lx,2π*n1/ly
            k = (kx^2 + ky^2)^0.5
            û[n1+ny,m1+1] = -im*ky/k^2*u[n1+ny,m1+1]
        end
    end
    û
end

xvelocityfield(d::AbstractDomain,u::Union{DSSField{T},GSSField{T}}) where {T<:AbstractFloat} = xvelocityfield(d,u.x[1])
xvelocity(d,u) = xvelocityfield(d,u) |> x->resolvedfield(d,x) |> x->inversefourier(d,x)

"""
    energyspectrum(d,u)
    Energy spectrum appropriately conjugated
"""
function scalarfield(d::AbstractDomain,u::DNSField{T};scalar='e') where {T<:AbstractFloat}
    (lx,ly),(nx,ny) = length(d),size(d)
    U = zeros(T,2ny-1,nx)
    for m1=0:nx-1
        nmin = m1==0 ? 1 : -ny+1
        for n1 = nmin:ny-1
            kx,ky = 2π*(m1/lx),2π*(n1/ly)
            k = scalar == 'e' ? (kx^2+ky^2)^0.5 : one(T)
            U[n1+ny,m1+1] = abs(u[n1+ny,m1+1])^2/k^2
        end
    end
    U
end

function scalarfield(d::AbstractDomain,u::DSSField{T};scalar='e') where {T<:AbstractFloat}
    (lx,ly),(nx,ny) = length(d),size(d)
    U = zeros(T,2ny-1,nx)
    m1 = 0
    for n1 = 1:ny-1
        kx,ky = 2π*(m1/lx),2π*(n1/ly)
        k = scalar == 'e' ? (kx^2+ky^2)^0.5 : one(T)
        U[n1+ny,m1+1] = abs(u.x[1][n1+ny])^2/k^2
    end
    for m1 = 1:nx-1
        for n1 = -ny+1:ny-1
            kx,ky = 2π*(m1/lx),2π*(n1/ly)
            k = scalar == 'e' ? (kx^2+ky^2)^0.5 : one(T)
            U[n1+ny,m1+1] =  abs(u.x[2][n1 + ny,n1 + ny,m1])/k^2
        end
    end
    U
end

function scalarfield(d::AbstractDomain,u::GSSField{T};scalar='e') where {T<:AbstractFloat}
    (lx,ly),(nx,ny),Λ = length(d),size(d),size(u.x[1],2)-1
    U = zeros(T,2ny-1,nx)
    for m1=0:Λ
        nmin = m1==0 ? 1 : -ny+1
        for n1 = nmin:ny-1
            kx,ky = 2π*(m1/lx),2π*(n1/ly)
            k = scalar == 'e' ? (kx^2+ky^2)^0.5 : one(T)
            U[n1+ny,m1+1] = abs(u.x[1][n1+ny,m1+1])^2/k^2
        end
    end
    for m1=Λ+1:nx-1
        for n1 = -ny+1:ny-1
            kx,ky = 2π*(m1/lx),2π*(n1/ly)
            k = scalar == 'e' ? (kx^2+ky^2)^0.5 : one(T)
            U[n1+ny,m1+1] = abs(u.x[2][n1+ny,m1-Λ,n1+ny,m1-Λ])/k^2
        end
    end
    U
end

energyspectrum(d,u) = scalarfield(d,u,scalar='e') |> x->resolvedfield(d,x)
enstrophyspectrum(d,u) = scalarfield(d,u,scalar='z') |> x->resolvedfield(d,x)
forcingspectrum(d,u0) = scalarfield(d,u0,scalar='e') |> x-> resolvedfield(d,x)

"""
    zonalenergy(d,u)
    zonalenstrophy(d,u)
    Zonal quadratic invariants for all modes including conjugates
"""
function zonalenergy(d::AbstractDomain,u)
    Ê = energyspectrum(d,u)
    E = [sum(Ê[:,d.nx])]
    append!(E,[2*sum(Ê[:,m]) for m=d.nx+1:2d.nx-1])
    E
end

function zonalenstrophy(d::AbstractDomain,u)
    Ẑ = enstrophyspectrum(d,u)
    Z = [sum(Ẑ[:,d.nx])]
    append!(Z,[2*sum(Ẑ[:,m]) for m=d.nx+1:2d.nx-1])
    Z
end

function modalenergy(d::AbstractDomain,u;m::Int=0)
    Ê = energyspectrum(d,u)
    Ê[:,d.nx+m]
end

"""
    energy(d,u)
    enstrophy(d,u)
    Quadratic invariants for all modes including conjugates
"""
energy(d::AbstractDomain,u) = sum(zonalenergy(d,u)) # the beauty of julia is killing me!
enstrophy(d::AbstractDomain,u) = sum(zonalenstrophy(d,u))

"""
    timeaverage(t,u)
    Time average anything given to it
"""
function timeaverage(t,u;t0=100.0)
    U = copy(u)
    t0 = min(t0,t[end])
    i0 = max(findfirst(x -> x > t0,t),2)
    for i=i0:length(u)
        U[i] = mean(u[i0-1:i])
    end
    U
end

"""
    modaleigvals(d,u)
    modalffts(d,u)
    Rank information for CE2
"""
function modaleigvals(d,u::DSSField{T}) where {T<:AbstractFloat}
    mEVs = zeros(T,2d.ny-1,d.nx-1)
    for m1=1:d.nx-1
        mEVs[:,m1] = real.(eigvals(u.x[2][:,:,m1]))
    end
    mEVs
end

secondcumulant(d,u::DSSField{T};m::Int=1) where {T<:AbstractFloat} = inversefourier(Domain(d.ly,d.ly,d.ny,d.ny),u.x[2][:,:,m])
secondcumulant(d,u;m::Int=1) = secondcumulant(d,convert(CE2(),u,d),m=m)

"""
    zonostrophy(d,u)
    Zonosotrophy and other invariants
"""
function zonostrophy(d,u::DNSField{T}) where {T<:AbstractFloat}
    E,Z = energy(length(d)...,size(d)...,u)
    U = sqrt(2E/(4π))
    ε = prob.c.μ*U^2
    β = 2prob.c.Ω*cos(prob.c.θ)
    LR = sqrt(2U/β)
    Lε = 0.5*(ε/β^3)^0.2
    LR/Lε
end

"""
    adjacency(d,u)
    Adjacency matrices for GQL/NL
"""
function adjacency(d::AbstractDomain;Λ=nx-1) where {T<:AbstractFloat}
    (lx,ly),(nx,ny) = length(prob.d),size(prob.d)
    B = bcoeffs(lx,ly,nx,ny,10.0,0.01,0.0,1.0) # set unity linear coefficients
    B̂ = zeros(Complex{T},2ny-1,2nx-1)
    for m=0:nx-1
        nmin = m==0 ? 1 : -ny+1
        for n = nmin:ny-1
            B̂[n+ny,m+nx] = B[n+ny,m+1]
            B̂[-n+ny,-m+nx] = conj(B[n+ny,m+1])
        end
    end

    M = 2nx-1
    N = 2ny-1
    Bij = zeros(Complex{T},M*N,M*N)

    for i=1:length(vec(B̂))
        Bij[i,i] = vec(B̂)[i]
    end

    Cij = zeros(Complex{T},M*N,M*N)
    Cp,Cm = ccoeffs(lx,ly,nx,ny,Λ)

    Ĉ = zeros(Complex{T},2ny-1,2nx-1,2ny-1,2nx-1)
    for m1=0:nx-1
        for n1=-ny+1:ny-1
            for m2=0:nx-1
                for n2=-ny+1:ny-1

                    Ĉ[n2+ny,m2+nx,n1+ny,m1+nx] = Cp[n2+ny,m2+1,n1+ny,m1+1]
                    Ĉ[-n2+ny,-m2+nx,-n1+ny,-m1+nx] = conj(Ĉ[n2+ny,m2+nx,n1+ny,m1+nx])

                    Ĉ[-n2+ny,-m2+nx,n1+ny,m1+nx] = Cm[n2+ny,m2+1,n1+ny,m1+1]
                    Ĉ[n2+ny,m2+nx,-n1+ny,-m1+nx] = conj(Ĉ[-n2+ny,-m2+nx,n1+ny,m1+nx])

                end
            end
        end
    end
    Cij = reshape(Ĉ,M*N,M*N)
    Bij,Cij
end

# container for solution based adjacency calculation
# to be changed to length u
function adjacency(d::AbstractDomain,u::Vector{DNSField{T}}) where {T<:AbstractFloat}
    U = zeros(Complex{T},(2d.nx-1)*(2d.ny-1),(2d.nx-1)*(2d.ny-1),1)
    for i=1:1
        @views U[:,:,i] .= adjacency(lx,ly,nx,ny,u[i])
    end
    U
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using JuMP
import JSON

Base.@kwdef mutable struct JuMPInstance <: Instance
    py::Union{Nothing,PyCall.PyObject} = nothing
    model::Union{Nothing,JuMP.Model} = nothing
    samples::Vector{PyCall.PyObject} = []

    function JuMPInstance(model::JuMP.Model)::JuMPInstance
        init_miplearn_ext(model)
        instance = new(nothing, model, [])
        py = PyJuMPInstance(instance)
        instance.py = py
        return instance
    end
end

function to_model(instance::JuMPInstance)::JuMP.Model
    return instance.model
end

function get_instance_features(instance::JuMPInstance)::Union{Vector{Float64},Nothing}
    return instance.model.ext[:miplearn]["instance_features"]
end

function _concat_features(dict, names)::Matrix{Float64}
    if isempty(dict)
        return zeros(length(names), 1)
    end
    ncols = length(first(dict).second)
    return vcat([n in keys(dict) ? dict[n]' : zeros(ncols) for n in names]...)
end

function _concat_categories(dict, names)::Vector{String}
    return String[n in keys(dict) ? dict[n] : n for n in names]
end

function get_variable_features(
    instance::JuMPInstance,
    names::Vector{String},
)::Matrix{Float64}
    return _concat_features(instance.model.ext[:miplearn]["variable_features"], names)
end

function get_variable_categories(instance::JuMPInstance, names::Vector{String})
    return _concat_categories(instance.model.ext[:miplearn]["variable_categories"], names)
end

function get_constraint_features(
    instance::JuMPInstance,
    names::Vector{String},
)::Matrix{Float64}
    return _concat_features(instance.model.ext[:miplearn]["constraint_features"], names)
end

function get_constraint_categories(instance::JuMPInstance, names::Vector{String})
    return _concat_categories(instance.model.ext[:miplearn]["constraint_categories"], names)
end

get_samples(instance::JuMPInstance) = instance.samples

function create_sample!(instance::JuMPInstance)
    sample = MemorySample()
    push!(instance.samples, sample)
    return sample
end

function find_violated_lazy_constraints(instance::JuMPInstance, solver)::Vector{String}
    if "lazy_find_cb" ∈ keys(instance.model.ext[:miplearn])
        return instance.model.ext[:miplearn]["lazy_find_cb"](instance.model, solver.data)
    else
        return []
    end
end

function enforce_lazy_constraint(instance::JuMPInstance, solver, violation::String)::Nothing
    instance.model.ext[:miplearn]["lazy_enforce_cb"](instance.model, solver.data, violation)
end

function solve!(solver::LearningSolver, model::JuMP.Model; kwargs...)
    solve!(solver, JuMPInstance(model); kwargs...)
end

function __init_PyJuMPInstance__()
    @pydef mutable struct Class <: miplearn.Instance
        function __init__(self, jl)
            self.jl = jl
        end
        to_model(self) = to_model(self.jl)
        get_instance_features(self) = get_instance_features(self.jl)
        get_variable_features(self, names) =
            get_variable_features(self.jl, from_str_array(names))
        get_variable_categories(self, names) =
            to_str_array(get_variable_categories(self.jl, from_str_array(names)))
        get_constraint_features(self, names) =
            get_constraint_features(self.jl, from_str_array(names))
        get_constraint_categories(self, names) =
            to_str_array(get_constraint_categories(self.jl, from_str_array(names)))
        get_samples(self) = get_samples(self.jl)
        create_sample(self) = create_sample!(self.jl)
        find_violated_lazy_constraints(self, solver, _) =
            find_violated_lazy_constraints(self.jl, solver)
        enforce_lazy_constraint(self, solver, _, violation) =
            enforce_lazy_constraint(self.jl, solver, violation)
    end
    copy!(PyJuMPInstance, Class)
end

export JuMPInstance, save, load_instance


text:
using HTsi
using Test
using Dates
using GMT
using RemoteFiles
using Printf

gmtrdt(dt) = Dates.format(dt, "yyyy-mm-ddTHH:MM:SS")

function plot_data(dts, data, fn="test")
    uts = datetime2unix.(dts)
    region = (gmtrdt(dts[1]), gmtrdt(dts[end]),
                minimum(data), maximum(data),)
    plot(uts, data,
        region=region,
        frame=(axes = :WSen,),
        xaxis=(annot = :auto, ticks = :auto),
        yaxis=(annot = :auto, ticks = :auto),
        fmt=:png,
        savefig=fn,
    )

end

function test_flux(enabled_files, dts)
    @testset "HTsi.jl fluxtable" begin
        force_download = false
        the_path = nothing

        the_path = path(HTsi._fluxtable)
        fluxtable = (enabled_files == nothing) || (:fluxtable in enabled_files)
        fluxtable && HTsi._init_fluxtable(local_path=the_path,
                                force_download=force_download)

        jds = datetime2julian.(dts)
        data = get_space_index.(F10obs(), jds)

        plot_data(dts, data, "test_flux")
    end
end

function test_dtcfile(enabled_files, dts)
    @testset "HTsi.jl dtcfile" begin

        force_download = false
        the_path = nothing

        # the_path = path(HTsi._dtcfile)
        dtcfile   = (enabled_files == nothing) || (:dtcfile in enabled_files)
        dtcfile && HTsi._init_dtcfile(local_path=the_path,
                            force_download=force_download)

        jds = datetime2julian.(dts)
        data = get_DstΔTc.(jds)

        plot_data(dts, data, "test_dtc")
    end
end

function test_solfsmy(enabled_files, dts)
    @testset "HTsi.jl solfsmy" begin

        force_download = false
        the_path = nothing

        # the_path = path(HTsi._solfsmy)
        solfsmy   = (enabled_files == nothing) || (:solfsmy in enabled_files)
        solfsmy && HTsi._init_solfsmy(local_path=the_path,
                            force_download=force_download)

        jds = datetime2julian.(dts)
        data = get_space_index.(S10(), jds)

        plot_data(dts, data, "test_solf")
    end
end

function test_wdcfiles(enabled_files, dts)
    @testset "HTsi.jl wdcfiles Kp Ap" begin

        force_download = false

        wdcfiles_oldest_year = 2000
        wdcfiles_newest_year = year(now()) - 1

        wdcfiles_dir = nothing

        wdcfiles  = (enabled_files == nothing) || (:wdcfiles in enabled_files)
        wdcfiles && HTsi._init_wdcfiles( force_download=force_download,
                            wdcfiles_oldest_year=wdcfiles_oldest_year,
                            wdcfiles_newest_year=wdcfiles_newest_year)

        # jds = datetime2julian.(dts)
        data = get_Kp.(dts)
        plot_data(dts, data, "test_wdc_kp")

    end
end

function test_dstfiles(enabled_files, dts)
    @testset "HTsi.jl dstfiles" begin
        force_download = false

        dstfiles_oldest_year = 2000
        dstfiles_newest_year = year(now()) - 1
        # dstfiles_newest_year = 2010
        # dts = DateTime(2010, ):Hour(1):DateTime(2010, 12)

        dstfiles_dir = nothing

        dstfiles  = (enabled_files == nothing) || (:wdcfiles in enabled_files)
        dstfiles && HTsi._init_dstfiles( force_download=force_download,
                            dstfiles_oldest_year=dstfiles_oldest_year,
                            dstfiles_newest_year=dstfiles_newest_year)

        # jds = datetime2julian.(dts)
        # data = get_space_index.(Dst(), jds)

        data = get_Dst.(dts)
        plot_data(dts, data, "test_dst")
    end
end

function test()
    dts = DateTime(2005, ):Day(1):DateTime(2019, )
    # dts = [DateTime(2015, 4, 25, 6, 11, 25),]
    enabled_files = nothing

    test_flux(enabled_files, dts)
    test_dtcfile(enabled_files, dts)
    test_solfsmy(enabled_files, dts)
    test_wdcfiles(enabled_files, dts)
    test_dstfiles(enabled_files, dts)

end
test()


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This code performs a convergence test for two finite difference methods applied to
Laplace's equation with variable diffusivity κ(x)."

# variable diffusivity
κ(x) = 1 + .5*sin(pi*x) # assumed κ(x) >= κ_min > 0
# κ(x) = Float64(1 + (x>0))

# method of manufactured solutions: define the exact solution, then figure out
# what f,α,β should be to make it the solution to your PDE
uexact(x) = log(2+sin(pi*x))
dudx_exact(x) = ForwardDiff.derivative(uexact,x)
f(x) = -ForwardDiff.derivative(x->κ(x)*dudx_exact(x),x) # -d/dx (κ(x)*du/dx) = f
α = uexact(-1)
β = uexact(1)

function solve(m,κ,f,α,β)
    x = LinRange(-1,1,m+2) # points x_0, x_1, ... x_m, x_m+1
    xint = x[2:end-1]
    h = x[2]-x[1]

    xmid = @. x[1:end-1] + h/2 # x_0 + h/2, x_1 + h/2, ..., x_m + h/2
    # (re-indexing using fractional indices) x_1/2, x_3/2, ..., x_i+1/2, ..., x_m+1/2
    k = κ.(xmid)

    #                        κ_{i-1/2}+κ_{i+1/2}      κ_{i+1/2}      κ_{i+1/2}
    A = (1/h^2) * spdiagm(0=>k[1:end-1]+k[2:end],-1=>-k[2:end-1],1=>-k[2:end-1])
    b = f.(xint)
    b[1] += k[1]*α/h^2
    b[end] += k[end]*β/h^2

    u = A\b
    return u,x,xint,h
end

m = 40
u,x,xint,h = solve(m,κ,f,α,β)
plot(xint,u,mark=:dot,markersize=2,label="FD solution")
plot!(x,uexact.(x),label="Exact solution")
plot!(leg=:bottomright)

# mvec = 2 .^ (2:9)
# hvec = zeros(length(mvec))
# err = zeros(length(mvec))
# for (i,m) in enumerate(mvec)
#     u,x,xint,h = solve(m,κ,f,α,β)
#     hvec[i] = h
#     err[i] = maximum(@. abs(uexact(xint) - u))
# end
# plot(hvec,err,marker=:dot,label="Max error")
# plot!(hvec,hvec.^2,linestyle=:dash,label="O(h^2)")
# plot!(xaxis=:log,yaxis=:log,legend=:bottomright)


text:
#=

Utilities for Continuous Data Augmentation in
Biogeographic Histories.

Ignacio Quintero Mächler

t(-_-t)

May 01 2017

=#

"""
    rejsam!(times::Array{Float64,1},
            si   ::Int64,
            sf   ::Int64,
            λ1   ::Float64,
            λ0   ::Float64,
            t    ::Float64)

Rejection sampling for an area conditioned on start and end point.
"""
function rejsam!(times::Array{Float64,1},
                 si   ::Int64,
                 sf   ::Int64,
                 λ1   ::Float64,
                 λ0   ::Float64,
                 t    ::Float64)

  samf = brprop!(times, si, λ1, λ0, t)

  while samf != sf
    samf = brprop!(times, si, λ1, λ0, t)
  end

  return nothing
end

"""
    brprop!(times::Array{Float64,1},
            si   ::Int64,
            λ1   ::Float64,
            λ0   ::Float64,
            t    ::Float64)

Two state DA proposal for one area.
"""
function brprop!(times::Array{Float64,1},
                 si   ::Int64,
                 λ1   ::Float64,
                 λ0   ::Float64,
                 t    ::Float64)
  empty!(times)

  c_st   ::Int64   = si
  c_time ::Float64 = 0.0
  endtime::Float64 = t

  re::Float64 = c_st == 0 ? rexp(λ1) : rexp(λ0)

  c_time += re

  while c_time < t
    push!(times, re)
    endtime  = t - c_time
    c_st     = 1 - c_st
    re       = c_st == 0 ? rexp(λ1) : rexp(λ0)
    c_time  += re
  end

  push!(times, endtime)

  return c_st
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA