text:
using LinearAlgebra
using SparseArrays
using DifferentialEquations
using Statistics

function matrix_multiplies(s)
    r = (0:1:10000)
    @time for n in 1:s
        r = test_matrix(3,.1)*r
    end
    return r
end

function lambda_func(t, p)
    base = p["sin_base"]
    rel_amp = p["sin_rel_amp"]
    rate = p["sin_rate"]
    phase_shift = p["phase_shift"]
    return base .+ base.*rel_amp.* cos.(rate*t.+phase_shift)
end

function bd_stat(Q)
    Q = Q'
    n = size(Q)[1]
    stat_dist = zeros(n)
    stat_dist[1]=1
    for i in range(2, stop=n)
        stat_dist[i] = stat_dist[i-1]*Q[i-1,i]/Q[i,i-1]
    end
    stat_dist = stat_dist./max(abs(sum(stat_dist)))
    if is_stochastic_vector(stat_dist)
        return stat_dist
    else
        print(stat_dist)
        print("NO!")
        return -1
    end
end

function bd_stat_mode(Q)
    Q = Q'
    n = size(Q)[1]
    births = diag(Q,1)
    deaths = diag(Q,-1)

    if size(findall(births./deaths.<1))[1]>0
        mode_index = minimum(findall(births./deaths.<1))
    else
        mode_index = n
    end
    stat_dist = zeros(n)
    stat_dist[mode_index] = 1

    for i in range(mode_index +1, stop = n)
        stat_dist[i] = stat_dist[i-1]*Q[i-1,i]/Q[i,i-1]
    end
    for i in reverse(range(1, stop=mode_index-1))
        stat_dist[i] = stat_dist[i+1]*Q[i+1,i]/Q[i,i+1]
    end
    stat_dist = stat_dist./max(abs(sum(stat_dist)))
    return stat_dist
    if is_stochastic_vector(stat_dist)
        return stat_dist
    else
        print("NO! \n")
        return -1
    end
end

function bd_stat_log(Q)
    Q = Q'
    n = size(Q)[1]
    stat_dist = zeros(n)
    stat_dist[1]=1
    stat_dist[1] = log(stat_dist[1])
    for i in range(2, stop=n)
        stat_dist[i] = stat_dist[i-1]+log.(Q[i-1,i])-log.(Q[i,i-1])
    end
    stat_dist = exp.(stat_dist)
    stat_dist = stat_dist./max(abs(sum(stat_dist)))
    if is_stochastic_vector(stat_dist)
        return stat_dist
    else
        return -1
    end
end

function bd_stat_mode_log(Q)
    Q = Q'
    n = size(Q)[1]
    births = diag(Q,1)
    deaths = diag(Q,-1)

    if size(findall(births./deaths.<1))[1]>0
        mode_index = minimum(findall(births./deaths.<1))
    else
        mode_index = n
    end
    stat_dist = zeros(n)
    stat_dist[mode_index] = 1
    stat_dist[mode_index] = log(stat_dist[mode_index])

    for i in range(mode_index +1, stop = n)
        stat_dist[i] = stat_dist[i-1]+log.(Q[i-1,i])-log.(Q[i,i-1])
    end
    for i in reverse(range(1, stop=mode_index-1))
        stat_dist[i] = stat_dist[i+1] + log.(Q[i+1,i])-log.(Q[i,i+1])
    end
    stat_dist = exp.(stat_dist)
    stat_dist = stat_dist./max(abs(sum(stat_dist)))
    return stat_dist
    if is_stochastic_vector(stat_dist)
        return stat_dist
    else
        print("NO! \n")
        return -1
    end
end

function Finite_Server_Q(t, p)
    S = Int(p["S"])
    s = Int(p["s"])
    mu = p["mu"]
    lambd = lambda_func(t,p)
    v = collect(1:S)
    br = ones(S)*lambd
    dr = min.(v, s)*mu
    dg = zeros(S+1)
    dg[1] = -br[1]
    dg[2:end-1] = -br[2:end]-dr[1:end-1]
    dg[end] = -dr[end]
    dr = float(dr)
    Q = Tridiagonal(dr,dg,br)
    return Q'
end

function Finite_Server_dudt(mu_t,p,t)
    if p["forwards"] == true
        return Finite_Server_Q(t, p)*mu_t
    else
        return Finite_Server_Q(1440.0-t, p)'*mu_t
    end
end

function Finite_Server_Number_Waiting_r(p)
    S = Int(p["S"])
    s = Int(p["s"])
    a = zeros(s)
    b = range(0, stop=S-s)
    return vcat(a,b)
end

function Infinite_Server_closed_form(S, t, lambd, mu)
    alpha = lambd/mu * (1-exp.(-mu*t))
    a = collect(1:S)
    a = pushfirst!(a,1)
    b = cumsum(log.(a))
    r = collect(0:1.0:S)
    mu_t = exp.(-alpha .+ log.(alpha)*r .- b)
    return mu_t
end

function Infinite_Server_forward(S, t, lambd, mu, mu_0)
    prob = ODEProblem(Infinite_Server_dudt, mu_0, (0.0,t),[mu, lambd, S])
    solve(prob, reltol=1e-8,abstol=1e-8,alg_hints = "stiff")    # alg_hints = "stiff",
end

function test_function(a,b)
    if a < 1 && b==3
        c = 5
    else
        c=0
    end
    return c
end

# US Callcenter Model

US_Bank_S = 470

function US_Bank_Lambda(t)
    if isa(t,Float64) || isa(t,Int)
        return US_Bank_Lambda_single(t)
    else
        t = collect(t)
        n = size(t)[1]
        lambdas = zeros(n)
        for i in range(1, stop=n)
            lambdas[i] = US_Bank_Lambda_single(t[i])
        end
        return lambdas
    end
end

function US_Bank_Lambda_single(t)

    if t <= 5.0*60.
      x_1, y_1 = 0.0,500.0/60.
      x_2, y_2 = 5.0*60.0,400.0/60.0

    elseif t > 5.0*60.0 && t <= 7.0*60.0
      x_1, y_1 = 5.0*60.0, 400.0/60.0
      x_2, y_2 = 7.0*60.0, 500.0/60.0

    elseif t > 7.0*60 && t <= 10.0*60
      x_1, y_1 = 7.0*60.0, 500.0/60
      x_2, y_2 = 10.0*60.0, 6000.0/60

    elseif t > 10.0*60 && t <= 16.0*60
      x_1, y_1 = 10.0*60., 6000.0/60
      x_2, y_2 = 16.0*60., 4500.0/60

    elseif t > 16.0*60 && t <= 18.0*60
      x_1, y_1 = 16.0*60., 4500.0/60
      x_2, y_2 = 18.0*60., 2000.0/60.

    elseif t > 18.0*60 && t <= 24.0*60
      x_1, y_1 = 18.0*60., 2000.0/60.
      x_2, y_2 = 24.0*60., 400.0/60.

    end
    m = (y_2 - y_1) / (x_2 - x_1)
    b = y_1 - m * x_1
    return m * t + b
end

function US_Bank_Servers_single(t)
    stdev = 3
    mu = 1.0/3.0
    if t <= 5. * 60.
        lbar = mean([500. / 60.,400. / 60])
    elseif t > 5.0 * 60 && t <= 7.0 * 60
        lbar = mean([400. / 60., 500.0 / 60.])
    elseif t > 7.0 * 60 && t <= 10.0 * 60
        lbar = mean([500.0 / 60, 6000.0 / 60])
    elseif t > 10.0 * 60 && t <= 16.0 * 60
        lbar = mean([6000.0 / 60.0, 4500.0 / 60.0])
    elseif t > 16.0 * 60 && t <= 18.0 * 60
        lbar = mean([4500.0 / 60., 2000. / 60.0])
    elseif t > 18.0 * 60 && t <= 24.0 * 60
        lbar =  mean([2000.0 / 60.0, 400.0 / 60.0])
    end
    return ceil(Int,lbar / mu + stdev * sqrt(lbar / mu))
end

function US_Bank_Servers(t)
    if isa(t,Float64) || isa(t,Int)
        return US_Bank_Servers_single(t)
    else
        t = collect(t)
        n = size(t)[1]
        servers = zeros(n)
        for i in range(1, stop=n)
            servers[i] = US_Bank_Servers_single(t[i])
        end
        return servers
    end
end

function US_Bank_Q(t)
    S = 470
    mu = 1/3.
    v = collect(1:S)
    br = ones(S)*US_Bank_Lambda(t)
    dr = min.(v, US_Bank_Servers(t)).*mu
    dg = zeros(S+1)
    dg[1] = -br[1]
    dg[2:end-1] = -br[2:end]-dr[1:end-1]
    dg[end] = -dr[end]
    dr = float(dr)
    Q = Tridiagonal(dr,dg,br)
    return Q'
end

function US_Bank_dudt(mu_t,p,t)
    if p["forwards"] == true
        return US_Bank_Q(t)*mu_t # mu, lambd
    else
        return US_Bank_Q(t)'*mu_t
    end
end

function US_Bank_forward()
    S = 470
    t = 1440
    mu_0 = zeros(S+1)
    mu_0[1] = 1
    prob = ODEProblem(US_Bank_dudt, mu_0, (0.0,t))
    solve(prob, reltol=1e-8, abstol=1e-8, alg_hints = "stiff")    # alg_hints = "stiff",
end

function is_Rate_Matrix(Q)
    if maximum(abs.(sum(Q, dims=2)))<1e-12 && all(diag(Q).<=0) && all(Q-Diagonal(diag(Q)).>=0)
        return true
    else
        return false
    end
end

function US_Bank_Forward_Euler(mu_0, t_f)
    h_sum = 0.0
    t = 0
    while(h_sum <= t_f)
        Qt = US_Bank_Q(t)
        h = 1/unif_param(Qt)
        mu_0 = mu_0 + (h/2.0)*Qt*mu_0
        h_sum += h
    end
    return mu_0, h_sum
end

function US_Bank_Number_Waiting_r(t)
    S = 470
    a = zeros(US_Bank_Servers(t))
    b = range(0, stop=S-US_Bank_Servers(t))
    return vcat(a,b)
end

function unif_param(Q)
    return maximum(abs.(diag(Q)))
end

function US_Bank_Backward_Euler(mu_0, t_f)
    h_sum = 0.0
    t = 0
    while(h_sum <= t_f)
        Qt = US_Bank_Q(t)
        h = 1/unif_param(Qt)
        mu_0 = mu_0 + (h/2.0)*Qt*mu_0
        h_sum += h
    end
    return mu_0, h_sum
end

function get_stat_dist_Q(Q)
    # ADD PETER WAY
    # if typeof(Q)== Tridiagonal{Float64,Array{Float64,1}}:
    #
    #
    # else
    n = size(Q)[1]
    B = Q[1:n-1, 1:n-1]
    d = Array{Float64,1}(Q[1:n-1,n])
    piv = B \ (-1*d)
    piv = vcat(piv,1)
    piv = piv/sum(piv)
    if ~ is_stochastic_vector(piv)
        piv[piv.<0] .= 0
        if is_stochastic_vector(piv)
            return piv
        else
            print("not stochastic after shift")
            return false
        end
    end
    return piv
end

function is_stochastic_matrix(P)
    truth = true
    if any(abs.(sum(P,dims=2).-1).>1e-7)
        truth = false
    elseif any(P.<0)
        truth = false
    elseif any(isnan.(P))
        truth = false
    elseif any(isinf.(P))
        truth = false
    end
    return truth
end

function is_stochastic_vector(piv)
    truth = true
    if abs(sum(piv)-1)>1e-8
        truth = false
    elseif ~all(piv.>=0)
        truth = false
    elseif any(isnan.(piv))
        truth = false
    elseif any(isinf.(piv))
        truth = false
    end
    return truth
end

function e_i(i,n)
    vec = zeros(n)
    vec[i]=1
    return vec
end

# Currently implemented for solution matrices
function tau(sol)
    return (maximum(sol,dims=1)-minimum(sol, dims=1))
end

function S_Selector(rho_max, eps=1e-8)
    S=50
    while((1-rho_max)*(rho_max^S)>eps)
        S = S+1
    end
    return S
end

function generate_primatives(chain_dict)
    if p["type"] == "sinusoidal"
        rho_max = chain_dict["rho_bar"] *  (1 + chain_dict["sin_rel_amp"])
        chain_dict["S"] = S_Selector(rho_max)
        chain_dict["sin_base"] = chain_dict["s"]*chain_dict["mu"]*chain_dict["rho_bar"]
        if chain_dict["sin_type"] == "rev-up"
            chain_dict["phase_shift"] = pi
            chain_dict["sin_rate"] = (1*pi)/1440.0
            print("dude")
        elseif chain_dict["sin_type"] == "Linda_Green"
            chain_dict["phase_shift"] = 0
            chain_dict["sin_rate"] = (2*pi)/1440.0
        elseif chain_dict["sin_type"] == "rev-up-four-times"
            chain_dict["phase_shift"] = pi
            chain_dict["sin_rate"] = (8*pi)/1440.0
        end
    end
    return chain_dict
end

# reward functions

function prob_of_wait_r(S, num_servers)
    a = zeros(num_servers+1)
    b = ones(S-num_servers)
    return vcat(a,b)
end

function number_waiting_r(S, num_servers)
    a = zeros(num_servers)
    b = range(0, stop=S-num_servers)
    return vcat(a,b)
end

max_stars_repo_name: alexinfanger/blogfiles

max_issues_repo_name: alexinfanger/blogfiles

max_forks_repo_name: alexinfanger/blogfiles


text:
using HMCExamples

HMCExamples.main_rbc_2_joint()

ext: jl


text:
function Base.convert( ::Type{Vector{String}},
                      x::DAQStringBuffer)
    length(x.str) == 0 && return nothing
    x.str[end] = 0
    strvec = split(unsafe_string(pointer(x.str)), ", "; keepempty=false)
    return String.(strvec)
end

function Base.convert(::Type{String}, x::DAQStringBuffer)
    x.str[end] = 0
    return unsafe_string(pointer(x.str))
end

function catch_error(code::Int32)
    code == DAQmx.Success && return
    buf = DAQStringBuffer(DAQmx.GetErrorString(code, Vector{UInt8}(),UInt32(0)))
    DAQmx.GetErrorString(code, buf.str, buf.size)
    code > 0 && (@warn convert(String, buf); return)
    code < 0 && throw(convert(String, buf))
end

function version()
    major, minor, update = Ref{UInt32}(), Ref{UInt32}(), Ref{UInt32}()
    DAQmx.GetSysNIDAQMajorVersion(major) |> catch_error
    DAQmx.GetSysNIDAQMinorVersion(minor) |> catch_error
    DAQmx.GetSysNIDAQUpdateVersion(update) |> catch_error
    major, minor, update = Int(major[]), Int(minor[]), Int(update[])
    return VersionNumber(major,minor,update)
end

function Base.show(io::IO, chan::PhysicalChannel{T}) where T <: AbstractIO
    println(io, "Channel name: ", chan.name)
    println(io, "Channel type: ", "$T")
    println(io, "Parent device: ", chan.parent.name)
    if chan.ranges !== nothing
        print(io, "Available ranges: ")
        for x in chan.ranges
            print(io, x, " ")
        end
    end
end

#FIXME: make a pretty print display for devices
function Base.show(io::IO, dev::DAQDevice)
    return nothing
end

function Base.push!(task::DAQTask{AI},
                      chan::PhysicalChannel{AI};
                     alias::String                  = "",
                      tcfg::DAQmx.DAQmxConstant     = DAQmx.Diff,
                     range::Tuple{Float64,Float64}  = (-10.0,10.0),
                     units::DAQmx.DAQmxConstant     = DAQmx.Volts,
                 scalename::String                  = "")
    #NOTE: TaskChannel not implemented yet--will simplify this...
    if (task.channels !== nothing && chan ∈ task.channels)
        return @warn "Channel $(chan.name) already appended to task."
    end
    DAQmx.CreateAIVoltageChan(task.handle, chan.name, alias,
                              tcfg, range[1], range[2],
                              units, scalename) |> catch_error
    task.device === nothing && (task.device = chan.parent)
    task.channels === nothing && (task.channels = [chan]; return)
    push!(task.channels, chan)
    return
end

function Base.push!(task::DAQTask{AI}, index::ChannelIndex, dev::DAQDevice=DefaultDev())
    if any(index .< 1)
        @error "Channel index not in valid range"
        return nothing
    end
    new_channels = dev.channels[AI][index]
    !isa(new_channels, Vector) && (push!(task, new_channels); return)
    for chan in new_channels
        push!(task, chan)
    end
end


text:
function lowpass(threshold::Float64)
    a(x) = x < threshold ? x : threshold
end

function highpass(threshold::Float64)
    a(x) = x > threshold ? x : threshold
end


text:
#=
Here we solve the equation:
```math
 q + dot(∇, uq) = 0
 p - dot(∇, up) = 0
```
on a sphere to test the conservation of the numerics

The boundary conditions are `p = q` when `dot(n, u) > 0` and
`q = p` when `dot(n, u) < 0` (i.e., `p` flows into `q` and vice-sersa).
=#

using MPI
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using Random

const _nstate = 2
const _q, _p = 1:_nstate
const stateid = (qid = _q, pid = _p)
const statenames = ("q", "p")

const _nauxstate = 3
# use a random velocity field
@inline function velocity_initialization!(vel, x, y, z)
  @inbounds begin
    vel[1] = rand()
    vel[2] = rand()
    vel[3] = rand()
  end
end

# physical flux function
@inline function flux!(F, Q, _, vel, _)
  @inbounds begin
    u, v, w = vel[1], vel[2], vel[3]
    q, p = Q[_q], Q[_p]

    F[1, _q], F[2, _q], F[3, _q] =  u*q,  v*q,  w*q
    F[1, _p], F[2, _p], F[3, _p] = -u*p, -v*p, -w*p
  end
end
# physical flux function
@inline function numerical_flux!(F, nM, QM, _, velM, QP, _, velP, _)
  @inbounds begin
    uM, vM, wM = velM[1], velM[2], velM[3]
    unM = nM[1] * uM + nM[2] * vM + nM[3] * wM
    uP, vP, wP = velP[1], velP[2], velP[3]
    unP = nM[1] * uP + nM[2] * vP + nM[3] * wP
    # since the velocity field is random we average here the two sides
    un = (unP + unM) / 2

    if un > 0
      F[_q], F[_p] = un * QM[_q], -un * QP[_p]
    else
      F[_q], F[_p] = un * QP[_q], -un * QM[_p]
    end
  end
end
@inline function numerical_boundary_flux!(F, nM, QM, _, vel, QP, _, _, _, _)
  @inbounds begin
    u, v, w = vel[1], vel[2], vel[3]
    un = nM[1] * vel[1] + nM[2] * vel[2] + nM[3] * vel[3]

    if un > 0
      F[_q], F[_p] = un * QM[_q], -un * QM[_q]
    else
      F[_q], F[_p] = un * QM[_p], -un * QM[_p]
    end
  end
end

function run(mpicomm, N, Nhorz, Rrange, timeend, DFloat, dt)
  ArrayType = Array

  topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DFloat,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = Topologies.cubedshellwarp
                                         )

  # spacedisc = data needed for evaluating the right-hand side function
  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = _nstate,
                           flux! = flux!,
                           numerical_flux! = numerical_flux!,
                           numerical_boundary_flux! = numerical_boundary_flux!,
                           auxiliary_state_length = _nauxstate,
                           auxiliary_state_initialization! =
                           velocity_initialization!,
                          )

  # This is a actual state/function that lives on the grid
  initialcondition!(Q, x, y, z, _...) = begin
    @inbounds Q[_q], Q[_p] = rand(), rand()
  end
  Q = MPIStateArray(spacedisc, initialcondition!)

  lsrk = LowStorageRungeKutta(spacedisc, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  sum0 = weightedsum(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e
  sum(Q₀) = %.16e""" eng0 sum0

  max_mass_loss = DFloat(0)
  max_mass_gain = DFloat(0)
  cbmass = GenericCallbacks.EveryXSimulationSteps(1) do
    cbsum = weightedsum(Q)
    max_mass_loss = max(max_mass_loss, sum0 - cbsum)
    max_mass_gain = max(max_mass_gain, cbsum - sum0)
    nothing
  end
  solve!(Q, lsrk; timeend=timeend, callbacks=(cbmass,))

  # Print some end of the simulation information
  engf = norm(Q)
  sumf = weightedsum(Q)
  @info @sprintf """Finished
  norm(Q)            = %.16e
  norm(Q) / norm(Q₀) = %.16e
  norm(Q) - norm(Q₀) = %.16e
  max mass loss      = %.16e
  max mass gain      = %.16e
  initial mass       = %.16e
  """ engf engf/eng0 engf-eng0 max_mass_loss max_mass_gain sum0
  max(max_mass_loss, max_mass_gain) / sum0
end

using Test
let
  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())
  mpicomm = MPI.COMM_WORLD
  if MPI.Comm_rank(mpicomm) == 0
    ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
    loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
    global_logger(ConsoleLogger(stderr, loglevel))
  else
    global_logger(NullLogger())
  end

  dt = 1e-4
  timeend = 100*dt

  polynomialorder = 4

  mpicomm = MPI.COMM_WORLD
  Nhorz = 4
  Rrange = 1.0:0.25:2.0

  for DFloat in (Float64,) #Float32)
    Random.seed!(0)
    delta_mass = run(mpicomm, polynomialorder, Nhorz, Rrange, timeend, DFloat, dt)
    @test abs(delta_mass) < 1e-15
  end
end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
module Day15
__revise_mode__ = :eval

using Advent: save_input

const my_input = [6,4,12,1,20,0,16]
const input_e1 = [0,3,6]

function simulate(input, rounds)
    last = input[end]
    times = Dict(x => i for (i, x) in enumerate(input))
    for curr in length(input):rounds - 1
        diff = curr - get!(times, last, curr)
        times[last] = curr
        last = diff
    end
    last
end

solve_a(input) = simulate(input, 2020)
solve_b(input) = simulate(input, 30000000)
solve(input = my_input) = solve_a(input), solve_b(input)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
include("./WilsoncloverFermion.jl")

abstract type Wilson_Dirac_operators{Dim} <: Dirac_operator{Dim} end

struct Wilson_Dirac_operator{Dim,T,fermion} <: Wilson_Dirac_operators{Dim}  where T <: AbstractGaugefields
    U::Array{T,1}
    boundarycondition::Vector{Int8}
    _temporary_fermi::Vector{fermion}
    γ::Array{ComplexF64,3}
    rplusγ::Array{ComplexF64,3}
    rminusγ::Array{ComplexF64,3}
    κ::Float64 #Hopping parameter
    r::Float64 #Wilson term
    hopp::Array{ComplexF64,1}
    hopm::Array{ComplexF64,1}
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    cloverterm::Union{Nothing,WilsonClover}
    verbose_print::Verbose_print
    _temporary_fermion_forCG::Vector{fermion}
    #verbose::Union{Verbose_1,Verbose_2,Verbose_3}
end

struct Wilson_Dirac_operator_evenodd{Dim,T,fermion} <: Dirac_operator{Dim}  where T <: AbstractGaugefields
    parent::Wilson_Dirac_operator{Dim,T,fermion}
    function Wilson_Dirac_operator_evenodd(D::Wilson_Dirac_operator{Dim,T,fermion}) where {Dim,T,fermion}
        new{Dim,T,fermion}(D)
    end
end

struct γ5D_Wilson_operator{Dim,T,fermion} <: γ5D_operator
    dirac::Wilson_Dirac_operator{Dim,T,fermion}
    function γ5D_Wilson_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
        return new{Dim,eltype(U),typeof(x)}(Wilson_Dirac_operator(U,x,parameters))
    end

    function γ5D_Wilson_operator(D::Wilson_Dirac_operator{Dim,T,fermion}) where {Dim,T,fermion}
        return new{Dim,T,fermion}(D)
    end
end

include("./WilsonFermion_4D.jl")

include("./WilsonFermion_4D_wing.jl")
include("./WilsonFermion_4D_nowing.jl")
include("./WilsonFermion_4D_wing_Adjoint.jl")

include("./WilsonFermion_2D.jl")
include("./WilsonFermion_2D_wing.jl")

#include("./WilsonFermion_4D_wing_fast.jl")

function Wilson_Dirac_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where {NC,Dim}
    xtype = typeof(x)
    num = check_parameters(parameters,"numtempvec",7)
    #num = 7
    _temporary_fermi = Array{xtype,1}(undef,num)

    @assert haskey(parameters,"κ") "parameters should have the keyword κ"
    κ = parameters["κ"]
    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end

    #boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])

    r = check_parameters(parameters,"r",1.0)

    if Dim==4
        γ,rplusγ,rminusγ = mk_gamma(r)
        hopp = zeros(ComplexF64,4)
        hopm = zeros(ComplexF64,4)
        hopp .= κ
        hopm .= κ
    elseif Dim == 2
        γ,rplusγ,rminusγ = mk_sigma(r)
        hopp = zeros(ComplexF64,2)
        hopm = zeros(ComplexF64,2)
        hopp .= κ
        hopm .= κ
    end

    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    numcg = check_parameters(parameters,"numtempvec_CG",7)
    #numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    #println("eps_CG = ",eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)
    verbose_print = Verbose_print(verbose_level)

    method_CG = check_parameters(parameters,"method_CG","bicg")

    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    hasclover = check_parameters(parameters,"hasclover",false)
    if hasclover
        cSW = r = check_parameters(parameters,"cSW",1.0)
        cloverterm = WilsonClover(cSW)
        #error("notsupported")
    else
        cloverterm = nothing
    end

    return Wilson_Dirac_operator{Dim,eltype(U),xtype}(U,boundarycondition,_temporary_fermi,
        γ,
        rplusγ,
        rminusγ,
        κ,
        r,
        hopp,
        hopm,
        eps_CG,MaxCGstep,verbose_level,
        method_CG,
        cloverterm,
        verbose_print,
        _temporary_fermion_forCG
        )
end

function (D::Wilson_Dirac_operator{Dim,T,fermion})(U) where {Dim,T,fermion}
    return Wilson_Dirac_operator{Dim,T,fermion}(U,D.boundarycondition,D._temporary_fermi,
        D.γ,
        D.rplusγ,
        D.rminusγ,
        D.κ,
        D.r,
        D.hopp,
        D.hopm,
        D.eps_CG,D.MaxCGstep,D.verbose_level,
        D.method_CG,
        D.cloverterm,
        D.verbose_print,
        D._temporary_fermion_forCG
        )
end

struct Adjoint_Wilson_operator{T} <: Adjoint_Dirac_operator
    parent::T
end

function Base.adjoint(A::T) where T <: Wilson_Dirac_operator
    Adjoint_Wilson_operator{typeof(A)}(A)
end

struct Adjoint_Wilson_operator_evenodd{T} <: Adjoint_Dirac_operator
    parent::T
end

function Base.adjoint(A::T) where T <: Wilson_Dirac_operator_evenodd
    Adjoint_Wilson_operator_evenodd{typeof(A)}(A)
end

function Initialize_WilsonFermion(u::AbstractGaugefields{NC,Dim};nowing = false) where {NC,Dim}
    _,_,NN... = size(u)
    return Initialize_WilsonFermion(NC,NN...,nowing = nowing)
end

function Initialize_4DWilsonFermion(u::AbstractGaugefields{NC,Dim}) where {NC,Dim}
    _,_,NN... = size(u)
    return WilsonFermion_4D_wing{NC}(NN...)
end

function Initialize_WilsonFermion(NC,NN...;nowing = false)
    Dim = length(NN)
    if Dim == 4
        if nowing
            fermion = WilsonFermion_4D_nowing{NC}(NN...)
        else
            fermion = WilsonFermion_4D_wing{NC}(NN...)
        end
        #fermion = WilsonFermion_4D_wing(NC,NN...)
    elseif Dim == 2
        fermion = WilsonFermion_2D_wing{NC}(NN...)
    else
        error("Dimension $Dim is not supported")
    end
    return fermion
end

using InteractiveUtils

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Wilson_Dirac_operator, T3 <:AbstractFermionfields}

    #@time Wx!(y,A.U,x,A,A._temporary_fermi)
    Wx!(y,A.U,x,A)
    if A.cloverterm != nothing
        Wclover!(y,A.U,x,A)
        error("not implemented!")
    end
    #error("w")
    #error("LinearAlgebra.mul!(y,A,x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")
end

function LinearAlgebra.mul!(y::T1,A::Wilson_Dirac_operator_evenodd{Dim,T,fermion},x::T3) where {T1 <:AbstractFermionfields,T, Dim,fermion, T3 <:AbstractFermionfields}
    WWx!(y,A.parent.U,x,A.parent)
end

function LinearAlgebra.mul!(y::T1,A::Adjoint_Wilson_operator_evenodd,x::T3) where {T1 <:AbstractFermionfields, T3 <:AbstractFermionfields}
    WWdagx!(y,A.parent.parent.U,x,A.parent.parent)
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_Wilson_operator, T3 <:  AbstractFermionfields}
    #error("LinearAlgebra.mul!(y,A,x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")
    Wdagx!(y,A.parent.U,x,A.parent)
    #error("LinearAlgebra.mul!(y,A,x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")

    return
end

"""
mk_gamma()
c----------------------------------------------------------------------c
c     Make gamma matrix
c----------------------------------------------------------------------c
C     THE CONVENTION OF THE GAMMA MATRIX HERE
C     ( EUCLIDEAN CHIRAL REPRESENTATION )
C
C               (       -i )              (       -1 )
C     GAMMA1 =  (     -i   )     GAMMA2 = (     +1   )
C               (   +i     )              (   +1     )
C               ( +i       )              ( -1       )
C
C               (     -i   )              (     -1   )
C     GAMMA3 =  (       +i )     GAMMA4 = (       -1 )
C               ( +i       )              ( -1       )
C               (   -i     )              (   -1     )
C
C               ( -1       )
C     GAMMA5 =  (   -1     )
C               (     +1   )
C               (       +1 )
C
C     ( GAMMA_MU, GAMMA_NU ) = 2*DEL_MU,NU   FOR MU,NU=1,2,3,4
c----------------------------------------------------------------------c
"""
function mk_gamma(r)
    g0 = zeros(ComplexF64,4,4)
    g1 = zero(g0)
    g2 = zero(g1)
    g3 = zero(g1)
    g4 = zero(g1)
    g5 = zero(g1)
    gamma = zeros(ComplexF64,4,4,5)
    rpg = zero(gamma)
    rmg = zero(gamma)

    g0[1,1]=1.0; g0[1,2]=0.0; g0[1,3]=0.0; g0[1,4]=0.0
    g0[2,1]=0.0; g0[2,2]=1.0; g0[2,3]=0.0; g0[2,4]=0.0
    g0[3,1]=0.0; g0[3,2]=0.0; g0[3,3]=1.0; g0[3,4]=0.0
    g0[4,1]=0.0; g0[4,2]=0.0; g0[4,3]=0.0; g0[4,4]=1.0

    g1[1,1]=0.0; g1[1,2]=0.0; g1[1,3]=0.0; g1[1,4]=-im
    g1[2,1]=0.0; g1[2,2]=0.0; g1[2,3]=-im;  g1[2,4]=0.0
    g1[3,1]=0.0; g1[3,2]=+im;  g1[3,3]=0.0; g1[3,4]=0.0
    g1[4,1]=+im;  g1[4,2]=0.0; g1[4,3]=0.0; g1[4,4]=0.0

    g2[1,1]=0.0; g2[1,2]=0.0; g2[1,3]=0.0; g2[1,4]=-1.0
    g2[2,1]=0.0; g2[2,2]=0.0; g2[2,3]=1.0; g2[2,4]=0.0
    g2[3,1]=0.0; g2[3,2]=1.0; g2[3,3]=0.0; g2[3,4]=0.0
    g2[4,1]=-1.0;g2[4,2]=0.0; g2[4,3]=0.0; g2[4,4]=0.0

    g3[1,1]=0.0; g3[1,2]=0.0; g3[1,3]=-im;  g3[1,4]=0.0
    g3[2,1]=0.0; g3[2,2]=0.0; g3[2,3]=0.0; g3[2,4]=+im
    g3[3,1]=+im;  g3[3,2]=0.0; g3[3,3]=0.0; g3[3,4]=0.0
    g3[4,1]=0.0; g3[4,2]=-im;  g3[4,3]=0.0; g3[4,4]=0.0

    g4[1,1]=0.0; g4[1,2]=0.0; g4[1,3]=-1.0;g4[1,4]=0.0
    g4[2,1]=0.0; g4[2,2]=0.0; g4[2,3]=0.0; g4[2,4]=-1.0
    g4[3,1]=-1.0;g4[3,2]=0.0; g4[3,3]=0.0; g4[3,4]=0.0
    g4[4,1]=0.0; g4[4,2]=-1.0;g4[4,3]=0.0; g4[4,4]=0.0

    g5[1,1]=-1.0;g5[1,2]=0.0; g5[1,3]=0.0; g5[1,4]=0.0
    g5[2,1]=0.0; g5[2,2]=-1.0;g5[2,3]=0.0; g5[2,4]=0.0
    g5[3,1]=0.0; g5[3,2]=0.0; g5[3,3]=1.0; g5[3,4]=0.0
    g5[4,1]=0.0; g5[4,2]=0.0; g5[4,3]=0.0; g5[4,4]=1.0

    gamma[:,:,1] = g1[:,:]
    gamma[:,:,2] = g2[:,:]
    gamma[:,:,3] = g3[:,:]
    gamma[:,:,4] = g4[:,:]
    gamma[:,:,5] = g5[:,:]

    for mu=1:4
        for j=1:4
            for i=1:4
                rpg[i,j,mu] = r*g0[i,j] + gamma[i,j,mu]
                rmg[i,j,mu] = r*g0[i,j] - gamma[i,j,mu]
            end
        end
    end

    return gamma,rpg,rmg

end

function mk_sigma(r)
    g0 = zeros(ComplexF64,2,2)
    g1 = zero(g0)
    g2 = zero(g1)
    g3 = zero(g1)

    gamma = zeros(ComplexF64,2,2,3)
    rpg = zero(gamma)
    rmg = zero(gamma)

    g0[1,1]=1.0; g0[1,2]=0.0
    g0[2,1]=0.0; g0[2,2]=1.0

    g1[1,1]=0.0; g1[1,2]=1.0
    g1[2,1]=1.0; g1[2,2]=0.0

    g2[1,1]=0.0; g2[1,2]=-im
    g2[2,1]=+im; g2[2,2]=0.0

    g3[1,1]=1.0; g3[1,2]=0.0
    g3[2,1]=0.0; g3[2,2]=-1.0

    gamma[:,:,1] = g1[:,:]
    gamma[:,:,2] = g2[:,:]
    gamma[:,:,3] = g3[:,:]

    for mu=1:2
        for j=1:2
            for i=1:2
                rpg[i,j,mu] = r*g0[i,j] + gamma[i,j,mu]
                rmg[i,j,mu] = r*g0[i,j] - gamma[i,j,mu]
            end
        end
    end

    return gamma,rpg,rmg

end

gtmp1,gtmp2,gtmp3 = mk_gamma(1)
const γ_all = gtmp1
const γ5 = γ_all[:,:,5]

const rplusγ1 = gtmp2
const rminusγ1 = gtmp3

include("./WilsontypeFermion.jl")

function Wx!(xout::T,U::Array{G,1},x::T,A,Dim)  where  {T,G <: AbstractGaugefields}
    #temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    #temp = temps[4]
    #temp1 = temps[1]cc
    #temp2 = temps[2]

    clear_fermion!(temp)
    #set_wing_fermion!(x)
    for ν=1:Dim

        xplus = shift_fermion(x,ν)
        #println(xplus)

        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication

        mul!(temp1,view(A.rminusγ,:,:,ν))

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)

        #
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rplusγ,:,:,ν),temp2)
        mul!(temp2,view(A.rplusγ,:,:,ν))

        add_fermion!(temp,A.hopp[ν],temp1,A.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)

    set_wing_fermion!(xout,A.boundarycondition)

    #display(xout)
    #    exit()
    return
end

function D4x!(xout::T1,U::Array{G,1},x::T2,A,Dim) where  {T1,T2,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    clear_fermion!(xout)
    set_wing_fermion!(x)

    for ν=1:Dim
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)
        #... Dirac multiplication
        mul!(temp1,view(A.rminusγ,:,:,ν),temp1)

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)
        mul!(temp2,Uminus',xminus)

        mul!(temp2,view(A.rplusγ,:,:,ν),temp2)
        add_fermion!(xout,0.5,temp1,0.5,temp2)

    end
    set_wing_fermion!(xout)
end

function D4dagx!(xout::T1,U::Array{G,1},x::T2,A,Dim) where  {T1,T2,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    clear_fermion!(xout)
    set_wing_fermion!(x)

    for ν=1:Dim
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)
        #... Dirac multiplication
        mul!(temp1,view(A.rplusγ,:,:,ν),temp1)

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)
        mul!(temp2,Uminus',xminus)

        mul!(temp2,view(A.rminusγ,:,:,ν),temp2)
        add_fermion!(xout,0.5,temp1,0.5,temp2)

    end
    set_wing_fermion!(xout)
end

function Dx!(xout::T1,U::Array{G,1},x::T2,A,Dim) where  {T1,T2,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    clear_fermion!(temp)
    #clear!(temp1)
    #clear!(temp2)
    set_wing_fermion!(x)
    for ν=1:Dim
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)
        #... Dirac multiplication
        mul!(temp1,view(A.rminusγ,:,:,ν),temp1)

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)
        mul!(temp2,Uminus',xminus)

        mul!(temp2,view(A.rplusγ,:,:,ν),temp2)
        add_fermion!(temp,0.5,temp1,0.5,temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1/(2*A.κ),x,-1,temp)

    #display(xout)
    #    exit()
    return
end

function Ddagx!(xout::T1,U::Array{G,1},x::T2,A,Dim) where  {T1,T2,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    clear_fermion!(temp)
    #clear!(temp1)
    #clear!(temp2)
    set_wing_fermion!(x)
    for ν=1:Dim
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)
        #... Dirac multiplication
        mul!(temp1,view(A.rplusγ,:,:,ν),temp1)

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)
        mul!(temp2,Uminus',xminus)

        mul!(temp2,view(A.rminusγ,:,:,ν),temp2)
        add_fermion!(temp,0.5,temp1,0.5,temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1/(2*A.κ),x,-1,temp)

    #display(xout)
    #    exit()
    return
end

function Tx!(xout::T,U::Array{G,1},x::T,A,Dim)  where  {T,G <: AbstractGaugefields} # Tx, W = (1 - T)x
    #temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]

    clear_fermion!(temp)
    #set_wing_fermion!(x)
    for ν=1:Dim

        xplus = shift_fermion(x,ν)
        #println(xplus)

        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication

        mul!(temp1,view(A.rminusγ,:,:,ν))

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)

        #
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rplusγ,:,:,ν),temp2)
        mul!(temp2,view(A.rplusγ,:,:,ν))

        add_fermion!(temp,A.hopp[ν],temp1,A.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,0,x,1,temp)

    set_wing_fermion!(xout,A.boundarycondition)

    #display(xout)
    #    exit()
    return
end

function Wdagx!(xout::T,U::Array{G,1},
    x::T,A,Dim) where  {T,G <: AbstractGaugefields}
    #,temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4] #temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    clear_fermion!(temp)
    #set_wing_fermion!(x)
    for ν=1:Dim
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication
        #mul!(temp1,view(x.rminusγ,:,:,ν),temp1)
        mul!(temp1,view(A.rplusγ,:,:,ν))

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rminusγ,:,:,ν),temp2)
        mul!(temp2,view(A.rminusγ,:,:,ν))

        add_fermion!(temp,A.hopp[ν],temp1,A.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)
    set_wing_fermion!(xout,A.boundarycondition)

    #display(xout)
    #    exit()
    return
end

include("./WilsonFermion_faster.jl")

struct DdagD_Wilson_operator{Dim,T,fermion} <: DdagD_operator
    dirac::Union{Wilson_Dirac_operator{Dim,T,fermion},Wilson_Dirac_operator_faster{Dim,T,fermion}}
    function DdagD_Wilson_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
        return new{Dim,eltype(U),typeof(x)}(Wilson_Dirac_operator(U,x,parameters))
    end

    function DdagD_Wilson_operator(D::Wilson_Dirac_operator{Dim,T,fermion}) where {Dim,T,fermion}
        return new{Dim,T,fermion}(D)
    end

    function DdagD_Wilson_operator(D::Wilson_Dirac_operator_faster{Dim,T,fermion}) where {Dim,T,fermion}
        return new{Dim,T,fermion}(D)
    end

end


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2

D = 64

#Q = Constant(randn(D,D))
#R = Constant(randn(D,D))
#Ψ = InfiniteCMPS(Q, R)

@load "DataQRD64_final" ΨL ρR HL E e hL

#alg = LBFGS(; verbosity = 2, maxiter = 10^6, gradtol = 1e-4);

Q = ΨL.Q
R = (ΨL.Rs)[1]

H = ∫(∂ψ[1]'*∂ψ[1] - 10000 * (ψ[1]'*ψ[1]) + 1000 * (ψ[1]'*ψ[1]'*ψ[1]*ψ[1]), (-Inf,+Inf))

ρL = one(ρR)
Z = tr(ρR)[]
_, E, e, hL, infoL = leftenv(H, (ΨL,ρL,ρR))
density = tr(R*ρR*R')[]

function environment(ΨL,ρL,ρR)
    linalg = GMRES(; tol = 1e-4)
    R = (ΨL.Rs)[1]
    density = tr(R*ρR*R')[]
    hL = R'*ρL*R
    hL = axpy!(-density, ρL, hL)

    HL₀ = zero(Q)

    let TL = LeftTransfer(ΨL)
        tol = linalg.tol
        HL, infoL = linsolve(hL, HL₀, linalg) do x
            y = ∂(x) - TL(x; tol = tol/10)
            y = axpy!(dot(ρR, x), ρL, y)
            CMPSKit.truncate!(y; tol = tol/10)
        end

        HL = rmul!(HL + HL', 0.5)
        HL = CMPSKit.truncate!(HL; tol = tol/10)
        return HL
    end
end
Env = environment(ΨL,ρL,ρR)
R² = dot(R,R)
@show Z
@show E
@show density
@show density² = density^2
@show tr(Env*ρR)
@show 2*tr(Env*R*ρR*R')[]
@show tr(R²*ρR*R²')
@show variance = sqrt(2*tr(Env*R*ρR*R')[] + density)
n̂ = ψ'*ψ
N̂ = ∫(n̂, (-Inf,+Inf));
@show N = expval(n̂, ΨL, ρL, ρR)
N = tr(R*ρR*R')[]
@show sqrt(2*tr(leftenv(N̂, (ΨL,ρL,ρR))[1]*R*ρR*R')[] + N)

@show psi = expval(ψ, ΨL, ρL, ρR)

# groundstate(H, Ψ; optalg = alg, linalg = GMRES(; tol = 1e-4), finalize! = finalize!)
#ΨL, ρR, E, e, normgrad, numfg, history = groundstate(H, Ψ; optalg = alg, linalg = GMRES(; tol = 1e-6))


text:
fname = :ST_ARCdense
c = Combi(hessian_dense,PDataST,solve_modelST_ARC,preprocessST,decreaseGen,TparamsST())
include("Template.jl")


text:
# Heta connector tests
heta_update()


text:
const BigInterval = IntervalArithmetic.Interval{BigFloat} #typealias BigInterval IntervalArithmetic.Interval{BigFloat}

function upbound(x)
  xu = _upbound(x)
  assert(xu≥0)
  xu
end
_upbound(x::Real) = x
_upbound(x::Interval) = x.hi

roundupbound(x) = ceil(Int,upbound(x))

prec{T}(::Type{Interval{T}}) = eps(T)
prec(T) = eps(T)

erroradjustment(T,err) = zero(T)
erroradjustment{T}(::Type{Interval{T}},err) = Interval(-upbound(err),upbound(err))
erroradjustment{IT<:Interval}(::Type{Complex{IT}},err) =
  complex(erroradjustment(IT,err),erroradjustment(IT,err))
erroradjustment(err) = erroradjustment(typeof(err),err)

intervalsign(S) =  (S > 0) ? 1 : (S < 0) ? -1 : error("$S contains zero")

# move to IntervalArithmetic
# import Base.log, Base.hypot
import Base.exp2, Base.sqrt

# VERSION < v"0.5" && (Base.hypot(x::Interval,y::Interval) = sqrt(x*x+y*y))
function Base.sqrt{T<:Interval}(z::Complex{T})
        (x,y) = reim(z)
        nz = hypot(x,y)
        if 0 ∈ y && x.lo < 0
          error("Not implemented")
          # if y == 0
          #   return 0∈x complex(y,sqrt(-x))
          # (xlh,ylh) = reim(csqrt(Interval(complex(x.lo,y.hi))))
          # (xhh,yhh) = reim(csqrt(Interval(complex(x.hi,y.hi))))
          # if y.lo == 0
          #   re = Interval(0,xhh.hi)
          #   im = Interval(sqrt(Interval(-x.hi)),xlh.hi)
          # else
          #   (xll,yll) = reim(csqrt(Interval(complex(x.lo,y.lo))))
          #   (xhl,yhl) = reim(csqrt(Interval(complex(x.hi,y.lo))))
          #   re = Interval(0,max(xhl.hi,xhh.hi))
          #   im = Interval(yll.lo,ylh.hi)
          # end
        else
          re = sqrt((nz+x)/2)
          im = y/2re
        end
        complex(re,im)
end
#
# Base.log{T<:Interval}(z::Complex{T}) = complex(log(abs(z)),atan2(imag(z),real(z)))
#
function Base.exp2{T<:Real}(z::Complex{T})
    er = exp2(real(z))
    theta = imag(z) * log(convert(T, 2))
    Complex(er*cos(theta), er*sin(theta))
end

# Akiyama–Tanigawa algorithm for second Bernoulli numbers B+n
function bernoulli2k(p::Integer,T)
  B2 = Array{T}(p+1)
  A = Array{BigInt}(2p+1)
  fact = one(BigInt)
  for m = 0:2p
      A[m+1] = fact
      fact *= m+1
      for j = m:(-1):1
          A[j] = j*((m+1)*A[j]-A[j+1])
      end
      iseven(m) && (B2[div(m,2)+1] = T(A[1])/fact
      )
  end
  return B2
end


text:
import Flux: Dense, softmax, σ, Chain

function foldlstyle()
    layers = [Dense(10,5,σ),Dense(5,2),softmax]
    model(x) = foldl((x,m) -> m(x), layers, init=x)
    @show model(rand(10))
end

function chain_style()
    model = Chain(Dense(10,5,σ),
                  Dense(5,2),
                  softmax)
    model(rand(10))
end

function circ_style()
    model = Dense(5,2)∘Dense(10,5,σ)
    model(rand(10))
end

function naive_function_chain()
    chain_func = Chain(x -> x^2, x -> x+1)
    chain_func(5)
end

@show foldlstyle()
@show chain_style()
@show circ_style()
@show naive_function_chain()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# indexing

export global_index, global_size, linear_index, @linearidx, @cartesianidx

## hardware

for f in (:blockidx, :blockdim, :threadidx, :griddim)
    @eval $f(ctx::AbstractKernelContext)::Int = error("Not implemented") # COV_EXCL_LINE
    @eval export $f
end

"""
    global_index(ctx::AbstractKernelContext)

Query the global index of the current thread in the launch configuration (i.e. as far as the
hardware is concerned).
"""
@inline function global_index(ctx::AbstractKernelContext)
    threadidx(ctx) + (blockidx(ctx) - 1) * blockdim(ctx)
end

"""
    global_size(ctx::AbstractKernelContext)

Query the global size of the launch configuration (total number of threads launched).
"""
@inline function global_size(ctx::AbstractKernelContext)
    griddim(ctx) * blockdim(ctx)
end

## logical

"""
    linear_index(ctx::AbstractKernelContext, grididx::Int=1)

Return a linear index for the current kernel by querying hardware registers (similar to
`get_global_id` in OpenCL). For applying a grid stride (in terms of [`global_size`](@ref)),
specify `grididx`.

"""
@inline function linear_index(ctx::AbstractKernelContext, grididx::Int=1)
    global_index(ctx) + (grididx - 1) * global_size(ctx)
end

"""
    linearidx(A, grididx=1, ctxsym=:ctx)

Macro form of [`linear_index`](@ref), which return from the surrouunding scope when out of
bounds:

    ```julia
    function kernel(ctx::AbstractKernelContext, A)
        idx = @linearidx A
        # from here on it's safe to index into A with idx
        @inbounds begin
            A[idx] = ...
        end
    end
    ```
"""
macro linearidx(A, grididx=1, ctxsym=:ctx)
    quote
        x = $(esc(A))
        i = linear_index($(esc(ctxsym)), $(esc(grididx)))
        i > length(x) && return
        i
    end
end

"""
    cartesianidx(A, grididx=1, ctxsym=:ctx)

Like [`@linearidx`](@ref), but returns a N-dimensional `CartesianIndex`.
"""
macro cartesianidx(A, grididx=1, ctxsym=:ctx)
    quote
        x = $(esc(A))
        i = @linearidx(x, $(esc(grididx)), $(esc(ctxsym)))
        @inbounds CartesianIndices(x)[i]
    end
end


text:
if !@isdefined(contains)
    contains(haystack, needle) = occursin(needle, haystack)
    export contains
end

# basic C stuff

puts(s::Union{String,SubString{String}}) = ccall(:puts, Cint, (Ptr{Cchar},), s)

jl_malloc(n::Integer) = ccall(:jl_malloc, Ptr{Cvoid}, (Csize_t,), n)

# check if a call failed

macro check(ex::Expr)
    ex.head == :call ||
        error("@check: not a call: $ex")
    arg1 = ex.args[1] :: Symbol
    if arg1 == :ccall
        arg2 = ex.args[2]
        arg2 isa QuoteNode ||
            error("@check: ccallee must be a symbol")
        f = arg2.value :: Symbol
    else
        f = arg1
    end
    prefix = "$f: "
    quote
        r = $(esc(ex))
        iszero(r) || @async @error($prefix * string(r))
        r
    end
end

# curl string list structure

struct curl_slist_t
    data::Ptr{Cchar}
    next::Ptr{curl_slist_t}
end


text:
dissolvedhandleref() = Reagents.Ref{Union{Catalyst,Nothing,Deleted}}(nothing)

struct DissolvedHandle
    ref::typeof(dissolvedhandleref())
end

DissolvedHandle() = DissolvedHandle(dissolvedhandleref())

function _tryset!(handle::DissolvedHandle, value)
    local old::Union{Catalyst{Any},Nothing,Deleted}  # `Catalyst{Any}` is OK for now
    old = handle.ref[]
    while true
        old isa Deleted && return old
        old, success = cas!(handle.ref, old, value)
        success && return old
    end
end

function Reagents.clear!(handle::DissolvedHandle)
    offer = _tryset!(handle, Deleted())
    offer isa Deleted && return
    offer === nothing && return
    _delete!(offer)
end

# TODO: make `DissolvedHandle` optional?
Reagents.dissolve(reagent::Reagent; once = false) =
    _dissolve(reagent, DissolvedHandle(), once)::DissolvedHandle

function _dissolve(reagent::Reagent, handle::DissolvedHandle, once::Bool = false)
    offer = Catalyst{Any}()  # TODO: narrow type
    let old = _tryset!(handle, offer)
        old isa Deleted && return handle
    end
    actr =
        then(once ? reagent : reagent ⨟ PostCommit(ReDissolve(reagent, handle)), Commit())
    ans = tryreact!(actr, nothing, Reaction(), offer)

    ntries = 0
    while true
        ans isa Block && return handle
        (once && !(ans isa Failure)) && return handle
        if ans isa NeedNack
            error("WithNack reagent cannot be dissolved")
            # ...or is it OK?
        end
        # Otherwise, it means that `reagent` just helped a reaction somewhere or
        # failed with `Retry`.  We need to make sure that there is no more duals
        # waiting for a swap.  However, since the catalyst is already dissolved,
        # we don't need to provide the `offer` here.
        ans = tryreact!(then(reagent, Commit()), nothing, Reaction(), nothing)

        ntries += 1
        should_limit_retries() && ntries > 1000 && error("too many retries")
    end
end
# An approach alternative to executing `redissolve` as above is to keep the
# catalyzing messages in the channel's bags. However, it'd require additional
# mechanism to re-execute pre-blocking reagents (e.g., `cas ⨟ swap`).

# Note: Registering catalysts on both ends of the channel can trigger an
# infinite loop. Is it possible to detect this?

struct ReDissolve{R<:Reagent}
    reagent::R
    handle::DissolvedHandle
end

(f::ReDissolve)(_) = _dissolve(f.reagent, f.handle)

redissolve!(_) = false
redissolve!(actr::Reactor) = redissolve!(actr.continuation)
function redissolve!(actr::Reactor{<:PostCommit{<:ReDissolve}})
    actr.reagent.f(nothing)
    return true
end

"""
    maybe_redissolve!(msg::Message)

This is called via `tryreact!(::Reactor{<:Swap}, ...)` when the reaction result
involving `msg` is `Retry()`. This could be due to that there is a stale CAS in
`msg.reaction`. We need to refresh these CASes to avoid deadlock.

TODO: better propagation of what CASes failed, to avoid reloading?
"""
function maybe_redissolve!(msg::Message)
    offer = msg.offer
    offer isa Catalyst || return
    has_stale_cas(msg.reaction) || return
    if cas!(offer.state, Pending(), Rescinded()).success
        redissolve!(msg.continuation)
    end
end


text:
type PlaneWaveMWTD{T,F,P} <: TDFunctional{T}
  direction::P
  polarisation::P
  speedoflight::T
  amplitude::F
end

# """
#     speedoflight(excitation)
#
# Returns the speed of light of the medium the excitation is defined in.
# """
# speedoflight(exc::PlaneWaveMWTD) = exc.speedoflight

function planewave(polarisation,direction,amplitude,speedoflight)
    PlaneWaveMWTD(direction,polarisation,speedoflight,amplitude)
end

*(a, pw::PlaneWaveMWTD) = PlaneWaveMWTD(
    pw.direction,
    a * pw.polarisation,
    pw.speedoflight,
    pw.amplitude
)

cross(k, pw::PlaneWaveMWTD) = PlaneWaveMWTD(
    pw.direction,
    k × pw.polarisation,
    pw.speedoflight,
    pw.amplitude
)

@compat function (f::PlaneWaveMWTD)(r,t)
    #dr = zero(typeof(t))
    t = cartesian(t)[1]
    #dr = zero(eltype(cartesian(r)))
    dr = zero(typeof(t))
    for i in 1 : 3
        dr += r[i]*f.direction[i]
    end
    f.polarisation * f.amplitude(t - dr/f.speedoflight)
end


text:
#
# These functions set up the structures needed for the MendelSearch routines,
# including the Parameter data structure whose entries control optimization.
#
export Parameter
export optimization_keywords!, set_parameter_defaults

mutable struct Parameter
  cases :: Int # number of cases in a least squares problem
  constraints :: Int # number of affine constraints
  goal :: AbstractString # "minimize" or "maximize"
  output_unit :: IO # for output of Search iterations
  parameters :: Int # number of parameters
  points :: Int # number of points in a grid search
  standard_errors :: Bool # true for parameter standard errors
  title :: AbstractString # problem title
  travel :: AbstractString # "grid" or "search"
  constraint :: Matrix{Float64} # parameter.constraint[constraint, par]
  constraint_level :: Vector{Float64}
  function_value :: Vector{Float64}
  grid :: Matrix{Float64} # parameter.grid[point, coordinate]
  min :: Vector{Float64}
  max :: Vector{Float64}
  name :: Vector{AbstractString}
  par :: Vector{Float64}
end

"""
Defines keywords specific to Search prior to parameter initialization.
"""
function optimization_keywords!(keyword::Dict{AbstractString, Any})

  keyword["cases"] = 0 # number of cases in a least squares problem
  keyword["constraints"] = 0 # number of affine constraints
  keyword["goal"] = "maximize" # "minimize" or "maximize"
  keyword["output_unit"] = stdout # for output of Search iterations
  keyword["parameters"] = 1 # number of parameters
  keyword["points"] = 0 # number of points in a grid search
  keyword["standard_errors"] = false # true for parameter standard errors
  keyword["title"] = "" # problem title
  keyword["travel"] = "search" # "search" or "grid"
  return keyword
end # optimization_keywords!

"""
Initializes the parameter data structure based on appropriate keywords.
"""
function set_parameter_defaults(keyword::Dict{AbstractString, Any})

  cases = keyword["cases"]
  constraints = keyword["constraints"]
  goal = keyword["goal"]
  output_unit = keyword["output_unit"]
  pars = keyword["parameters"]
  points = keyword["points"]
  standard_errors = keyword["standard_errors"]
  title = keyword["title"]
  travel = keyword["travel"]
  #
  # Check that input variables are properly defined.
  #
  if constraints < 0
    throw(ArgumentError(
      "The number of constraints ($constraints) must be nonnegative.\n \n"))
  end
  if goal != "minimize" && goal != "maximize"
    throw(ArgumentError(
      "The variable goal ($goal) must be minimize or maximize.\n \n"))
  end
  if pars < 0
    throw(ArgumentError(
      "The number of parameters ($pars) must be nonnegative.\n \n"))
  end
  if points < 0
    throw(ArgumentError(
      "The number of grid points ($points) must be nonnegative.\n \n"))
  end
  if travel != "grid" && travel != "search"
    throw(ArgumentError(
      "The variable travel ($travel) must be grid or search.\n \n"))
  end
  #
  # Name the parameters par 1, par 2, and so forth.
  #
  pname = ["par" for i = 1:pars]
  for i = 1:pars
    pname[i] = pname[i]*" $i"
    pname[i] = rpad(pname[i], 8, ' ')
  end
  #
  # Set other array defaults.
  #
  grid = zeros(points, pars)
  par = zeros(pars)
  pmin = zeros(pars)
  pmin[1:pars] .= -Inf
  pmax = zeros(pars)
  pmax[1:pars] .= Inf
  constraint = zeros(constraints, pars)
  constraint_level = zeros(constraints)
  level = zeros(constraints)
  if travel == "grid"
    function_value = zeros(points)
  else
    function_value = zeros(2)
  end
  parameter = Parameter(cases, constraints, goal, output_unit, pars,
    points, standard_errors, title, travel, constraint, constraint_level,
    function_value, grid, pmin, pmax, pname, par)
  #
  # These initial values can now be altered via the returned parameter.
  #
  return parameter
end # function set_parameter_defaults

max_stars_repo_name: OpenMendel/Search

max_issues_repo_name: OpenMendel/Search

max_forks_repo_name: OpenMendel/Search


text:
# ----------------------
# Holds simulations
# ----------------------

"Container holds simulations of all sections"
struct SharedSimulations{R<:AbstractRange}
    @addStructFields(SharedMatrix{Float64}, d0, d1, d0psi, d1psi, d0eur, d1eur, d0eursq, d1eursq, d0eurcub, d1eurcub)
    @addStructFields(SharedMatrix{Float64}, epsdeq1, epsdgt1, Prdeq1, Prdgt1)
    @addStructFields(SharedMatrix{Float64}, Eeps, profit, surplus, revenue, drillcost, extension)
    @addStructFields(SharedMatrix{Float64}, D_at_T)
    zchars_time::R
end

"Create `SharedSimulations`"
function SharedSimulations(pids::Vector{<:Integer}, nT::Integer, N::Integer, Dmax::Integer, zchars_time::AbstractRange)

    @declareVariables( SharedMatrix{Float64}(nT, N, pids=pids), d0, d1, d0psi, d1psi, d0eur, d1eur, d0eursq, d1eursq, d0eurcub, d1eurcub)
    @declareVariables( SharedMatrix{Float64}(nT, N, pids=pids), epsdeq1, epsdgt1, Prdeq1, Prdgt1)
    @declareVariables( SharedMatrix{Float64}(nT, N, pids=pids), Eeps, profit, surplus, revenue, drillcost, extension)

    D_at_T = SharedMatrix{Float64}(Dmax, N, pids=pids)

    return SharedSimulations(
        d0, d1, d0psi, d1psi, d0eur, d1eur, d0eursq, d1eursq, d0eurcub, d1eurcub,
        epsdeq1, epsdgt1, Prdeq1, Prdgt1,
        Eeps, profit, surplus, revenue, drillcost, extension,
        D_at_T,
        zchars_time
    )
end

function SharedSimulations(pids, data::DataDrill)
    nT = length(zchars(data))
    N = num_i(data)
    Dmax = _Dmax(statespace(_model(data)))
    zchars_time = _timestamp(zchars(data))
    return SharedSimulations(pids, nT, N, Dmax, zchars_time)
end

SharedSimulations(data::DataDrill) = SharedSimulations(workers(), data)

function fill!(s::SharedSimulations, x)
    fill!(s.d0        , x)
    fill!(s.d1        , x)
    fill!(s.d0psi     , x)
    fill!(s.d1psi     , x)
    fill!(s.d0eur     , x)
    fill!(s.d1eur     , x)
    fill!(s.d0eursq   , x)
    fill!(s.d1eursq   , x)
    fill!(s.d0eurcub  , x)
    fill!(s.d1eurcub  , x)
    fill!(s.epsdeq1   , x)
    fill!(s.epsdgt1   , x)
    fill!(s.Prdeq1    , x)
    fill!(s.Prdgt1    , x)
    fill!(s.Eeps      , x)
    fill!(s.profit    , x)
    fill!(s.surplus   , x)
    fill!(s.revenue   , x)
    fill!(s.drillcost , x)
    fill!(s.extension , x)
    fill!(s.D_at_T    , x)
    return nothing
end


text:
module TestCategories
using Test

using Catlab.Theories
using Catlab.CategoricalAlgebra.Sets, Catlab.CategoricalAlgebra.Categories

# Categories from Julia types
#############################

C = TypeCat(FreeCategory.Ob, FreeCategory.Hom)
@test Ob(C) == TypeSet(FreeCategory.Ob)
@test sprint(show, C) == "TypeCat($(FreeCategory.Ob), $(FreeCategory.Hom))"

F = id(C)
@test (dom(F), codom(F)) == (C, C)
@test startswith(sprint(show, F), "id(TypeCat(")

x, y = Ob(FreeCategory, :x, :y)
f = Hom(:f, x, y)
@test F(x) == x
@test F(f) == f

α = id(F)
@test (dom(α), codom(α)) == (F, F)
@test component(α, x) == id(F(x))

end


text:
export Nested

struct Nested{T}
    value::T
end

# function Nested(nt)
#     x = keysort(nt)
#     T = typeof(x)
#     Nested{T}(x)
# end

Nested(;kwargs...) = Nested((;kwargs...))

import Base

function Base.show(io::IO, x::Nested)
    io = IOContext(io, :compact => true)
    print(io, "Nested(", x.value, ")")
end

lenses(x::Nested) = lenses(x.value)

function schema(N::Type{Nested{X}}) where {X}
    return Nested(schema(X))
end

@generated function Base.empty(x::Nested)
    e = _empty(schema(x))
    return e
end

function _empty(x::Nested)
    ℓ = lenses(x)
    n = length(ℓ)
    set(x, batch(ℓ...), ntuple(i -> □, n))
end


text:
module Wigner

using StaticArrays
import ACE.Rotations3D.rotation_D_matrix

# D_{μm}^l(Ang2Mat_zyz(α,β,γ))
function wigner_big_D(μ,m,l,α,β,γ)
	return (exp(-im*α*m) * wigner_d(m,μ,l,β)  * exp(-im*γ*μ))'
end

# Wigner small d, modified from
# https://github.com/cortner/SlaterKoster.jl/blob/
# 8dceecb073709e6448a7a219ed9d3a010fa06724/src/code_generation.jl#L73
function wigner_d(μ, m, l, β)
    fc1 = factorial(l+m)
    fc2 = factorial(l-m)
    fc3 = factorial(l+μ)
    fc4 = factorial(l-μ)
    fcm1 = sqrt(fc1 * fc2 * fc3 * fc4)

    cosb = cos(β / 2.0)
    sinb = sin(β / 2.0)

    p = m - μ
    low  = max(0,p)
    high = min(l+m,l-μ)

    temp = 0.0
    for s = low:high
       fc5 = factorial(s)
       fc6 = factorial(l+m-s)
       fc7 = factorial(l-μ-s)
       fc8 = factorial(s-p)
       fcm2 = fc5 * fc6 * fc7 * fc8
       pow1 = 2 * l - 2 * s + p
       pow2 = 2 * s - p
       temp += (-1)^(s+p) * cosb^pow1 * sinb^pow2 / fcm2
    end
    temp *= fcm1

    return temp
end

function mat2ang(Q)
	return mod(atan(Q[2,3],Q[1,3]),2pi), acos(Q[3,3]), mod(atan(Q[3,2],-Q[3,1]),2pi);
end

# Rotation D matrix
function rot_D(L::Integer, Q)
	Mat_D = zeros(ComplexF64, 2L + 1, 2L + 1);
	D = rotation_D_matrix(L);
	α, β, γ = mat2ang(Q);
	for i = 1 : 2L + 1
		for j = 1 : 2L + 1
			Mat_D[i,j] = wigner_big_D(D[i,j].μ, D[i,j].m, L, α, β, γ);
		end
	end
	return SMatrix{2L+1, 2L+1, ComplexF64}(Mat_D)
end

function rand_QD(L)
	rotz(α) = [cos(α) -sin(α) 0; sin(α) cos(α) 0; 0 0 1]
	roty(α) = [cos(α) 0 sin(α); 0 1 0;-sin(α) 0 cos(α)]
	Ang2Mat_zyz(α,β,γ) = rotz(α)*roty(β)*rotz(γ)

	α = 2pi*rand();
	β = pi*rand();
	γ = 2pi*rand();

	# construct the Q matrix
	Q = Ang2Mat_zyz(α,β,γ)
	Q = SMatrix{3,3}(Q)

	rand_ref = rand((-1,1))

	return rand_ref * Q, rand_ref^L * rot_D(L, Q)
end

function rand_QD(L1, L2)
	rotz(α) = [cos(α) -sin(α) 0; sin(α) cos(α) 0; 0 0 1]
	roty(α) = [cos(α) 0 sin(α); 0 1 0;-sin(α) 0 cos(α)]
	Ang2Mat_zyz(α,β,γ) = rotz(α)*roty(β)*rotz(γ)

	α = 2pi*rand();
	β = pi*rand();
	γ = 2pi*rand();

	# construct the Q matrix
	Q = Ang2Mat_zyz(α,β,γ)
	Q = SMatrix{3,3}(Q)

	rand_ref = rand((-1,1))

	return rand_ref * Q, rand_ref^(L1) * rot_D(L1, Q), rand_ref^(L2) * rot_D(L2, Q)
end

end


text:
############# data preparation ###################

function parsetitle(title, r)
    run = r.data
    d = Dict{Symbol, Any}(k => missing for k in (:displace_location, :displace_direction , :nest_coverage, :transfer))
    d[:nest_coverage] = "open"
    d[:nest] = run.originalnest
    d[:feeder] = run.feeder
    d[:fictive_nest] = run.nest
    d[:track] = run.track
    d[:title] = title
    d[:comment] = r.metadata.comment
    for kv in split(title, ' ')
        k, v = split(kv, '#')
        if k ∉ ("person", "pellet")
            if k == "nest"
                d[:nest_coverage] = v
            else
                d[Symbol(k)] = v
            end
        end
    end
    return (; pairs(d)...)
end

function getdf(data)
    data["nest#closed person#therese displace_direction#none displace_location#feeder"] = data["nest#closed person#therese"]
    delete!(data, "nest#closed person#therese")
    df = DataFrame(parsetitle(k, r) for (k, v) in data for r in v.runs)

    @. df[!, :displace_direction] = switchdirections.(df.displace_direction)
    @. df[!, :group] .= getgroup(df.displace_location, df.transfer, df.displace_direction)
    @. df[!, :set] = getset(df.transfer, df.group)

    categorical!(df, [:group, :set, :displace_direction, :displace_location, :nest_coverage, :transfer])
    levels!(df.group, ["none", "left", "right", "away", "towards", "zero", "back", "far"])

    filter!(r -> r.group ≠ "far" || r.title == "transfer#far person#therese", df)

    df[!, :direction_deviation]  = [angle(r.fictive_nest - r.feeder, turningpoint(r.track) - r.feeder) for r in eachrow(df)]
    max_direction_deviation = maximum(r.direction_deviation for r in eachrow(df) if r.group ∉ ("far", "zero"))
    mean_direction_deviation = mean(r.direction_deviation for r in eachrow(df) if r.group ∉ ("far", "zero"))
    filter!(r -> r.group ≠ "far" || r.direction_deviation < 4mean_direction_deviation, df)

    df[!, :turning_point] .= zero.(df.feeder)
    df[!, :center_of_search] .= zero.(df.feeder)
    for r in eachrow(df)
        trans = createtrans(r.nest, r.displace_location, r.fictive_nest, r.feeder)
        @. r.track.coords = trans(r.track.coords)
        @. r.track.rawcoords.xy .= trans(r.track.rawcoords.xy)
        r.feeder = trans(r.feeder)
        r.fictive_nest = trans(r.fictive_nest)
        r.nest = trans(r.nest)
        Δ = intended(r.group) - r.fictive_nest
        r.turning_point = turningpoint(r.track) + Δ
        r.center_of_search = searchcenter(r.track) + Δ
    end

    groups = levels(df.group)
    nc = length(groups)
    colors = OrderedDict(zip(groups, [colorant"black"; distinguishable_colors(nc - 1, [colorant"white", colorant"black"], dropseed = true)]))

    gdf = groupby(df, :group)
    DataFrames.transform!(gdf, :group => (g -> colors[g[1]]) => :groupcolor)
    DataFrames.transform!(gdf, :groupcolor => getcolor => :color)

    df
end

switchdirections(_::Missing) = missing
switchdirections(d) =   d == "left" ? "right" :
                        d == "right" ? "left" :
                        d

getgroup(displace_location::Missing, transfer, displace_direction) = transfer
getgroup(displace_location, transfer, displace_direction) = displace_location == "nest" ? "zero" : displace_direction
getset(_::Missing, d) = d == "none" ? "Closed" : "Displacement"
getset(_, __) = "Transfer"

intended(d::AbstractString) =  d == "none" ? DungBase.Point(0,0) :
                d == "away" ? DungBase.Point(0, -50) :
                d == "towards" ? DungBase.Point(0, 50) :
                d == "right" ? DungBase.Point(50, 0) :
                d == "left" ? DungBase.Point(-50, 0) :
                d == "zero" ? DungBase.Point(0, -130) :
                d == "back" ? DungBase.Point(0,0) :
                d == "far" ? DungBase.Point(0,0) :
                error("unknown displacement")
intended(d) = intended(string(d))

_get_rotation_center(displace_location::Missing, nest, fictive_nest) = fictive_nest
_get_rotation_center(displace_location, nest, fictive_nest) = displace_location == "feeder" ? fictive_nest : nest
_get_zeroing(nest::Missing, fictive_nest) = fictive_nest
_get_zeroing(nest, fictive_nest) = nest
function createtrans(nest, displace_location, fictive_nest, feeder)
    v = feeder - _get_rotation_center(displace_location, nest, fictive_nest)
    α = atan(v[2], v[1])
    rot = LinearMap(Angle2d(-π/2 - α))
    trans = Translation(-_get_zeroing(nest, fictive_nest))
    passmissing(rot ∘ trans)
end

function highlight(c, i, n)
    h = HSL(c)
    HSL(h.h, h.s, i/(n + 1))
end
function getcolor(g)
    n = length(g)
    [highlight(c, i, n) for (i, c) in enumerate(g)]
end


text:
using BinaryBuilder

# These are the platforms built inside the wizard
platforms = [
    BinaryProvider.Linux(:i686, :glibc),
  BinaryProvider.Linux(:x86_64, :glibc),
  BinaryProvider.Linux(:aarch64, :glibc),
  BinaryProvider.Linux(:armv7l, :glibc),
  BinaryProvider.Linux(:powerpc64le, :glibc),
  BinaryProvider.MacOS(),
  BinaryProvider.Windows(:i686),
  BinaryProvider.Windows(:x86_64)
]

# If the user passed in a platform (or a few, comma-separated) on the
# command-line, use that instead of our default platforms
if length(ARGS) > 0
    platforms = platform_key.(split(ARGS[1], ","))
end
info("Building for $(join(triplet.(platforms), ", "))")

# Collection of sources required to build cfitsio
sources = [
    "http://heasarc.gsfc.nasa.gov/FTP/software/fitsio/c/cfitsio3420.tar.gz" =>
    "6c10aa636118fa12d9a5e2e66f22c6436fb358da2af6dbf7e133c142e2ac16b8",
]

script = raw"""
cd $WORKSPACE/srcdir
curl https://heasarc.gsfc.nasa.gov/docs/software/fitsio/cfileio.c > cfitsio/cfileio.c
cd cfitsio/
./configure --prefix=/ --host=$target --enable-reentrant
make -j10 shared
make -e install

"""

products = prefix -> [
    LibraryProduct(prefix,"libcfitsio")
]

# Build the given platforms using the given sources
hashes = autobuild(pwd(), "cfitsio", platforms, sources, script, products)

if !isempty(get(ENV,"TRAVIS_TAG",""))
    print_buildjl(pwd(), products, hashes,
        "https://github.com/Keno/CFITSIOBuilder/releases/download/$(ENV["TRAVIS_TAG"])")
end

max_stars_repo_name: Keno/CFITSIOBuilder

max_issues_repo_name: Keno/CFITSIOBuilder

max_forks_repo_name: Keno/CFITSIOBuilder


text:
export formula_FR, formula_PR, formula_HS, formula_HZ

function formula_FR(∇f, ∇ft, s, d)

    β = (∇ft⋅∇ft) / (∇f⋅∇f)

    return β
end

function formula_PR(∇f, ∇ft, s, d)

    y = ∇ft - ∇f
    β = (∇ft⋅y) / (∇f⋅∇f)

    return β
end

function formula_HS(∇f, ∇ft, s, d)

    y = ∇ft - ∇f
    β = (∇ft⋅y) / (d⋅y)

    return β
end

function formula_HZ(∇f, ∇ft, s, d)

    y = ∇ft - ∇f
    n2y = y⋅y
    β1 = (y⋅d)
    β = ((y - 2 * d * n2y / β1)⋅∇ft) / β1

    return β
end


text:
using PyCall
pygui(:qt)
using PyPlot
matplotlib[:rcdefaults]()
ion()

using Unwrap

x = linspace(-1,1,200)
y = x
r² = broadcast(+, x.^2, (y').^2)
z = exp(-r²)
φ = angle(exp(im*z*10π))

φ_old = copy(φ)
φ1 = unwrap2d(φ)
φ2 = unwrap2d(φ1) # Iterating improves result

function plot_map(args...;kwargs...)
    p = pcolormesh(args...;
                   rasterized=true,
                   cmap=plt[:cm][:get_cmap]("viridis"),
                   kwargs...)
    margins(0,0)
    p
end

function plot_phase(φ,gs,i,label,ylabels=false)
    fig[:add_subplot](get(gs, (1,i-1)))
    pm = plot_map(x,y,φ/2π)
    xlabel(L"x")
    if ylabels
        ylabel(L"y")
    else
        gca()[:set_yticklabels]([])
    end
    cax=fig[:add_subplot](get(gs, (0,i-1)))
    colorbar(pm,cax=cax,orientation="horizontal",label=label)
    cax[:set_xticklabels](cax[:get_xticklabels](),rotation=45)
    cax[:xaxis][:tick_top]()
    cax[:xaxis][:set_ticks_position]("both")
    cax[:xaxis][:set_label_position]("top")
end
fig = figure("2d",figsize=(10,4.5))
clf()
gs = matplotlib[:gridspec][:GridSpec](2,3,height_ratios=(1,10))
plot_phase(φ_old,gs,1,L"Wrapped [$\pi$]",true)
plot_phase(φ1,gs,2,L"1st unwrapped [$\pi$]")
plot_phase(φ2,gs,3,L"2nd unwrapped [$\pi$]")
gs[:tight_layout](fig)

savefig("2d.png",dpi=72,transparent=true)


text:
module Exercism

using JSON

"Open the notebook for the excersice 'slug' in the current directory and create a submittable 'slug.jl' file containing the user solution."
function create_submission(slug)
    open("$slug.jl", "w") do f
        write(f, parse_notebook("$slug.ipynb"))
    end
end

"Parses the notebook in 'path' and returns the code of all cells marked with # submit"
function parse_notebook(path)
    nb = open(JSON.parse, path, "r")

    # check if notebook is in acceptable format
    nb["nbformat"] == 4 || error("unrecognized notebook format ", nb["nbformat"])
    lang = lowercase(nb["metadata"]["language_info"]["name"])
    lang == "julia" || error("unrecognized notebook language $lang")

    submission = ""

    # scan for cells that contain code marked for submission with `# submit`
    for cell in nb["cells"]
        if cell["cell_type"] == "code" && !isempty(cell["source"])
            s = join(cell["source"])
            if startswith(s, "# submit")
                submission *= s
            end
        end
    end

    replace(submission, "# submit\n" => "")
end

end # module


text:
"""
# ElectricalEngineering.Gray

Defines colors based on the four base types of gray `colorBlack1`,
`colorBlack2`, `colorBlack3` and `colorBlack4`. So independent of the picked
color everything will always appear in gray scale. This way a gray scale and a
color plot can alternatively be created. So either `ElectricalEngineering.Gray` or `ElectricalEngineering.Tab20bc` can
be loaded (at a time) and be used.

# Colors

The following colors are substitued by `colorBlack1`,
`colorBlack2`, `colorBlack3` and `colorBlack4`:

- `colorBrown1`, `colorBrown2`, `colorBrown3`, `colorBrown4`
- `colorBlue1`, `colorBlue2`, `colorBlue3`, `colorBlue4`
- `colorRed1`, `colorRed2`, `colorRed3`, `colorRed4`
- `colorOrange1`, `colorOrange2`, `colorOrange3`, `colorOrange4`
- `colorGreen1`, `colorGreen2`, `colorGreen3`, `colorGreen4`
- `colorPurple1`, `colorPurple2`, `colorPurple3`, `colorPurple4`
- `colorDarkPurple1`, `colorDarkPurple2`, `colorDarkPurple3`, `colorDarkPurple4`
- `colorGray1`, `colorGray2`, `colorGray3`, `colorGray4`
- `colorMagenta1`, `colorMagenta2`, `colorMagenta3`, `colorMagenta4`

# Example

Copy and paste code:

```julia
using PyPlot,ElectricalEngineering,ElectricalEngineering.Gray
x=collect(0:0.01:1);
figure(figsize=(3.3,2.5))
plot(x,x,color=colorRed1,linestyle=lineStyle1,linewidth=lineWidth1,
label="(1)")
plot(x,x.^2,color=colorRed2,linestyle=lineStyle2,linewidth=lineWidth2,
label="(2)")
plot(x,x.^3,color=colorRed3,linestyle=lineStyle3,linewidth=lineWidth3,
label="(3)")
plot(x,x.^4,color=colorRed4,linestyle=lineStyle4,linewidth=lineWidth4,
label="(4)")
subplots_adjust(left=0.20, right=0.95, top=0.95, bottom=0.20)
legend(loc="upper left",fontsize=legendFontSize)
```

# See also

`ElectricalEngineering.Tab20bc`

"""
module Gray
    using ElectricalEngineering
    export colorBrown1,colorBrown2,colorBrown3,colorBrown4,
        colorBlue1,colorBlue2,colorBlue3,colorBlue4,
        colorRed1,colorRed2,colorRed3,colorRed4,
        colorOrange1,colorOrange2,colorOrange3,colorOrange4,
        colorGreen1,colorGreen2,colorGreen3,colorGreen4,
        colorPurple1,colorPurple2,colorPurple3,colorPurple4,
        colorDarkPurple1,colorDarkPurple2,colorDarkPurple3,colorDarkPurple4,
        colorMagenta1,colorMagenta2,colorMagenta3,colorMagenta4,
        colorGray1,colorGray2,colorGray3,colorGray4

    # Color scheme is only based on colorBlack1..4
    const colorBrown1=colorBlack1
    const colorBrown2=colorBlack2
    const colorBrown3=colorBlack3
    const colorBrown4=colorBlack4

    const colorBlue1=colorBlack1
    const colorBlue2=colorBlack2
    const colorBlue3=colorBlack3
    const colorBlue4=colorBlack4

    const colorRed1=colorBlack1
    const colorRed2=colorBlack2
    const colorRed3=colorBlack3
    const colorRed4=colorBlack4

    const colorOrange1=colorBlack1
    const colorOrange2=colorBlack2
    const colorOrange3=colorBlack3
    const colorOrange4=colorBlack4

    const colorGreen1=colorBlack1
    const colorGreen2=colorBlack2
    const colorGreen3=colorBlack3
    const colorGreen4=colorBlack4

    const colorPurple1=colorBlack1
    const colorPurple2=colorBlack2
    const colorPurple3=colorBlack3
    const colorPurple4=colorBlack4

    const colorDarkPurple1=colorBlack1
    const colorDarkPurple2=colorBlack2
    const colorDarkPurple3=colorBlack3
    const colorDarkPurple4=colorBlack4

    const colorMagenta1=colorBlack1
    const colorMagenta2=colorBlack2
    const colorMagenta3=colorBlack3
    const colorMagenta4=colorBlack4

    const colorGray1=colorBlack1
    const colorGray2=colorBlack2
    const colorGray3=colorBlack3
    const colorGray4=colorBlack4
end


text:
#################################################################################################
#                               Load Packages
#################################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("..")
using Revise, PAAV, Plots, DataFrames, Statistics, StatsPlots
using GLM, Random
include("simulation ratio effect.jl")
Random.seed!(524184)
#################################################################################################
#                               Conjunctive Search
#################################################################################################
set_sizes = [1,2,5,10,15]

conj_results = run_simulation(set_sizes; fun=conjunctive_ratio,
    Δτ=.39, topdown_weight=.66)

pyplot()
@df conj_results plot(:distractors, :hit_rate, grid=false,
    ylims=(.5,1), leg=false,ylabel="Accuracy", xlabel="Set Size",
    color=:grey, linewidth=2, xaxis=font(10), yaxis=font(10), size=(600,400))

#savefig("Conjunctive_hit_rate.eps")

@df conj_results plot(:distractors, :rt_mean, group=(:target_present,:response), grid=false,
     ylims=(0,3.0), linewidth=2, leg=true, ylabel="Mean RT (seconds)", xlabel="Set Size",
     size=(600,400))

#savefig("Conjunctive_set_size.eps")

df_present = filter(x->x[:target_present] ==:present && x[:response] ==:present, conj_results)
conj_ols_present = lm(@formula(rt_mean ~ distractors), df_present)

df_absent = filter(x->x[:target_present] ==:absent && x[:response] ==:absent, conj_results)
conj_ols_absent = lm(@formula(rt_mean ~ distractors), df_absent)

max_stars_repo_path: validation/Ratio Effect Validation.jl

max_issues_repo_path: validation/Ratio Effect Validation.jl

max_forks_repo_path: validation/Ratio Effect Validation.jl


text:
"""
    Return the graph whose vertices are the edges of `g`, which are linked by an edge if they share a vertex.
"""
function line_graph(g::AbstractSimpleGraph)
    E = incidence_matrix(g)
    A = E'*E - 2I
    return SimpleGraph(A)
end


text:
using Pkg
Pkg.activate(@__DIR__)

using ProgressiveHedging
import JuMP
import Ipopt

function two_stage_model(scenario_id::ScenarioID)

    model = JuMP.Model(()->Ipopt.Optimizer())
    JuMP.set_optimizer_attribute(model, "print_level", 0)
    JuMP.set_optimizer_attribute(model, "tol", 1e-12)
    JuMP.set_optimizer_attribute(model, "acceptable_tol", 1e-12)

    scen = value(scenario_id)

    ref = JuMP.@variable(model, x >= 0.0)
    stage1 = [ref]

    ref = JuMP.@variable(model, y >= 0.0)
    stage2 = [ref]

    b_s = scen == 0 ? 11.0 : 4.0
    c_s = scen == 0 ? 0.5 : 10.0

    JuMP.@constraint(model, x + y == b_s)

    JuMP.@objective(model, Min, 1.0*x + c_s*y)

    return JuMPSubproblem(model,
                          scenario_id,
                          Dict(stid(1) => stage1,
                               stid(2) => stage2)
                          )
end

scen_tree = two_stage_tree(2)

function my_callback(ext::Dict{Symbol,Any},
                     phd::PHData,
                     winf::ProgressiveHedging.WorkerInf,
                     niter::Int)
    # The `ext` dictionary can be used to store things between PH iterations
    if niter == 2
        ext[:message] = "This is from iteration 2!"
    elseif niter == 5
        println("Iteration 5 found the message: " * ext[:message])
    elseif niter == 10
        println("This is iteration 10!")
        # We can access the current consensus variable values
        for (xhid, xhat) in pairs(consensus_variables(phd))
            println("The value of $(name(phd,xhid)) is $(value(xhat)).")
        end
    end
    # Returning false from the callback will terminate PH.
    # Here we stop after 20 iterations.
    return niter < 20
end

my_cb = Callback(my_callback)

(niter, abs_res, rel_res, obj, soln_df, phd) = solve(scen_tree,
                                                     two_stage_model,
                                                     ScalarPenaltyParameter(1.0),
                                                     callbacks=[my_cb]
                                                     )
@show niter
@show abs_res
@show rel_res
@show obj
@show soln_df


text:
# Autogenerated wrapper script for Usrsctp_jll for x86_64-apple-darwin
export libusrsctp

JLLWrappers.@generate_wrapper_header("Usrsctp")
JLLWrappers.@declare_library_product(libusrsctp, "@rpath/libusrsctp.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libusrsctp,
        "lib/libusrsctp.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function fan_columns!(y::Matrix{Float64}, x::Matrix{Float64}, columns::Vector{Int64},
                      offset_x::Int64)
    my, ny = size(y)
    mx, nx = size(x)
    source = offset_x*mx + 1
    @inbounds for i = 1:length(columns)
        k = columns[i]
        destination = (k - 1)*my + 1
        copyto!(y, destination, x, source, mx)
        source += mx
    end
end

function get_abc!(a::Matrix{Float64}, b::Matrix, c::Matrix{Float64}, jacobian::Matrix{Float64}, m::Model)
    fill!(a, 0.0)
    fill!(b, 0.0)
    fill!(c, 0.0)
    i_bkwrd = m.i_bkwrd_b
    i_current = m.i_current
    i_fwrd = m.i_fwrd_b
    fan_columns!(a, jacobian, m.i_bkwrd_b, 0)
    offset = m.n_bkwrd + m.n_both
    fan_columns!(b, jacobian, m.i_current, offset)
    offset += m.n_current
    fan_columns!(c, jacobian, m.i_fwrd_b, offset)
end

"""
    function h0!(h0::Matrix{Float64}, a::Matrix{Float64}, b::Matrix{Float64}, c::Matrix{Float64})

computes h0 = inv(a+b*c)
"""
function h0!(h0::AbstractMatrix{Float64}, a::AbstractMatrix{Float64},
             b::AbstractMatrix{Float64}, c::AbstractMatrix{Float64},
             work::AbstractMatrix{Float64}, linsolve_ws)
    @inbounds copy!(work, a)
    @inbounds mul!(work, b, c, 1.0, 1.0)
    fill!(h0, 0.0)
    n = size(h0, 1)
    m = 1
    @inbounds for i = 1:n
        h0[m] = 1.0
        m += n + 1
    end
    @inbounds linsolve_core!(work, h0, linsolve_ws)
    if any(isnan.(h0))
        if any(isnan.(a))
            throw(ArgumentError("NaN in a"))
        end
        if any(isnan.(b))
            throw(ArgumentError("NaN in b"))
        end
        if any(isnan.(c))
            throw(ArgumentError("NaN in c"))
        end
        throw(ArgumentError("NaN in h0"))
    end
end

"""
function hh!(hh::AbstractMatrix{Float64}, h::AbstractMatrix{Float64},
             f::AbstractMatrix{Float64}, hf::AbstractMatrix{Float64},
             n::Int64, work1::AbstractMatrix{Float64},
             work2::AbstractMatrix{Float64})

computes hh = [h, h(1), h(2), ..., h(n-1)] and h(i) = -h*f*h(-1)
"""
function hh!(hh::AbstractMatrix{Float64}, h::AbstractMatrix{Float64},
             f::AbstractMatrix{Float64}, hf::AbstractMatrix{Float64},
             preconditioner_window::Int64, work1::AbstractMatrix{Float64},
             work2::AbstractMatrix{Float64})
    if any(isnan.(h))
        throw(ArgumentError("NaN in h"))
    end
    if any(isnan.(f))
        throw(ArgumentError("NaN in f"))
    end
    m = size(h, 1)
    m2 = m*m
    @inbounds mul!(hf, h, f, -1, 0)
    @inbounds copy!(work1, h)
    k = 1
    @inbounds for i = 1:preconditioner_window
        copyto!(hh, k, work1, 1, m2)
        if i < preconditioner_window
            mul!(work2, hf, work1)
            copy!(work1, work2)
            k += m2
        end
    end
end


text:
@everywhere using StanSample, Distributed

ProjDir = @__DIR__

themodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

# Default for tmpdir is to create a new tmpdir location
# To prevent recompilation of a Stan progam, choose a fixed location,
# and to initally hav access to run logs etc.
tmpdir = "$(ProjDir)/tmp"

N = 8
@time all_models = pmap(i -> SampleModel("sm_$i", themodel,
    tmpdir=tmpdir,
    method=StanSample.Sample(save_warmup=true, num_warmup=1000,
      num_samples=1000, adapt=StanSample.Adapt(delta=0.85))
  ), 1:N)

# This could be any Vector type (or Dict with Int keys)
observed_data = Dict(:N => 10, :y => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# This might need to be encapsulated (or contain all data for all models?)
@time rc = pmap(m->stan_sample(m; data=observed_data), all_models);

p = []
if all(success.(rc))
  for m in all_models
    push!(p, read_samples(m, output_format=:particles))
  end
  display(p)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetIdentity <: SwaggerModel
    principalId::Any # spec type: Union{ Nothing, String } # spec name: principalId
    tenantId::Any # spec type: Union{ Nothing, String } # spec name: tenantId
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function VirtualMachineScaleSetIdentity(;principalId=nothing, tenantId=nothing, type=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetIdentity, Symbol("principalId"), principalId)
        setfield!(o, Symbol("principalId"), principalId)
        validate_property(VirtualMachineScaleSetIdentity, Symbol("tenantId"), tenantId)
        setfield!(o, Symbol("tenantId"), tenantId)
        validate_property(VirtualMachineScaleSetIdentity, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type VirtualMachineScaleSetIdentity

const _property_map_VirtualMachineScaleSetIdentity = Dict{Symbol,Symbol}(Symbol("principalId")=>Symbol("principalId"), Symbol("tenantId")=>Symbol("tenantId"), Symbol("type")=>Symbol("type"))
const _property_types_VirtualMachineScaleSetIdentity = Dict{Symbol,String}(Symbol("principalId")=>"String", Symbol("tenantId")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ VirtualMachineScaleSetIdentity }) = collect(keys(_property_map_VirtualMachineScaleSetIdentity))
Swagger.property_type(::Type{ VirtualMachineScaleSetIdentity }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetIdentity[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetIdentity }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetIdentity[property_name]

const _allowed_VirtualMachineScaleSetIdentity_type = ["SystemAssigned"]

function check_required(o::VirtualMachineScaleSetIdentity)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetIdentity }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "VirtualMachineScaleSetIdentity", :enum, val, _allowed_VirtualMachineScaleSetIdentity_type)
    end
end


text:
module TestRotate

using Test

import TreeLas: Cycle
import GraphIdx
import GraphIdx: IncidenceIndex, Vec
import GraphIdx.Tree: ChildrenIndex, hierarchy_string
import GraphIdx.Tree: lowest_common_ancestors, root_node
import GraphIdx.Tree: dfs_finish

@testset "Rotate Grid(2, 3)              " begin
    g = GraphIdx.Grid.GridGraph(2, 3)
    if false
        GraphIdx.enumerate_edges(g) do i, u, v, _
            println(i, ":", u, " ", v)
        end
    end

    pi = [1, 1, 1, 3, 3, 5]
    cidx = ChildrenIndex(pi)
    @test hierarchy_string(cidx) == """
1
├─2
└─3
  ├─4
  └─5
    └─6
"""
    n = GraphIdx.num_nodes(g)
    m = GraphIdx.num_edges(g)
    lam = Vec(Float64.(collect(1:m)))
    cb = Cycle.CycleBasis(g, pi, cidx)
    @test cb.non_tree_edges == [2 => 4, 4 => 6]
    @test cb.lca == [1, 3]

    tlam0, tlam = Cycle.extract_rotate(cb, lam)
    @test tlam[2:end] ≈ [7, 10, 15, 12, 10]
    @test tlam0[2:end] ≈ [1, 4, 2, 5, 3]
end

@testset "Rotate Grid(2, 3) 2            " begin
    g = GraphIdx.Grid.GridGraph(2, 3)
    pi = [1, 1, 1, 6, 3, 5]
    cidx = ChildrenIndex(pi)
    @test hierarchy_string(cidx) == """
1
├─2
└─3
  └─5
    └─6
      └─4
"""
    n = GraphIdx.num_nodes(g)
    m = GraphIdx.num_edges(g)
    lam = Vec(Float64.(collect(1:m)))
    cb = Cycle.CycleBasis(g, pi)
    @test cb.non_tree_edges == [3 => 4, 2 => 4]
    @test cb.lca == [3, 1]

    tlam0, tlam = Cycle.extract_rotate(cb, lam)
    @test tlam0[2:end] ≈ [1, 4, 7, 5, 3]
    @test tlam[2:end] ≈ (tlam0[2:end] .+ 6) + [0, 0, 2, 2, 2]
end

@testset "Rotate Demo 3x7                " begin
    include("demo3x7.jl")
    n1, n2 = size(y)
    graph = GraphIdx.Grid.GridGraph(n1, n2)
    @test size(graph) == (n1, n2)
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
module ANOVADDPTest

using DataFrames
using Distributions
using DPMNeal3
using ExtractMacro
using Random
using SpecialFunctions
using StatsBase
using StatsModels

import Distributions: pdf, logpdf, shape, scale, rate, rand
import DPMNeal3: parent_dpm, logpredlik, update_hyperpars!, update_suffstats!

export BernoulliDDP, NormalDDP, PoissonDDP
export BernoulliData, NormalData, PoissonData
export train, simple_effect_probabilities
export interaction_effect_probabilities, expandgrid
export gamma_codebook, gamma_posterior, anova_bnp_normal
export anova_bnp_poisson
export group_codes, group_probs, effects1, effects2, fpost

include("gamma_posterior.jl")
include("gamma_codebook.jl")
include("expandgrid.jl")
include("womack.jl")
include("normalinversegamma.jl")
include("normal.jl")
include("poisson.jl")
include("bernoulli.jl")
include("dpm_mcmc_chain.jl")
include("train.jl")
include("simple_effect_probabilities.jl")
include("interaction_effect_probabilities.jl")
include("anova_bnp_normal.jl")
include("anova_bnp_poisson.jl")
include("anova_bnp_bernoulli.jl")

end # module


text:
module Tabben

export TabularDataset

include("dataset.jl")

end


text:
# import Dyno
#
# sys_str = Dict{String,String}(
#     "a"=>"1.0",
#     "b"=>"a+1",
#     "c"=>"a+b",
# )
#
# # res = Dyno.string_dict_to_symbol_dict(sys_str)
# res = Dyno.triangular_system(sys_str)
# #
# # sys = Dict(
# #     :a=>:(1.0),
# #     :b=>:(a+1),
# #     :c=>:(a+b),
# # )


text:
# DT can be block sparse structure or normal tensor
struct MPS{DT}
	Mvec::Vector{DT}

end

import Base.iterate
iterate(mps::MPS) = iterate(mps.Mvec)
iterate(mps::MPS, state) = iterate(mps.Mvec, state)

import Base.getindex
getindex(mps, idx::Int) = mps.Mvec[idx]

import Base.setindex
function setindex(mps::MPS, v, idx::Int)
	mps.Mvec[idx] = v
end

sysSize(mps::MPS) = length(mps.Mvec)
#
#	simple arithmetic
#

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS