text:
"""
    body_thrust(Ts, orientations; rotors_on = true)

Returns thrust in the body frame in 3-D.
"""
@inline function body_thrust(Ts, orientations; rotors_on = true)
    Ts_vec = Vector{ST.SArray{Tuple{3},Float64,1,3}}(undef,length(Ts)) # 2-D
    # Ts_vec = Vector{ST.SArray{Tuple{3},Float64,1,3}}(undef,length(Ts)) # 3-D
    if rotors_on
        for (i,T) in enumerate(Ts)
            # T_vec = ST.@SVector [orientations[i][1] * T, orientations[i][3]] # 2-D
            Ts_vec[i] = ST.@SVector [orientations[i][1] * T, orientations[i][2] * T, orientations[i][3] * T] # 3-D
        end
    else
        for (i,T) in enumerate(Ts)
            Ts_vec[i] = ST.@SVector [0.0, 0.0, 0.0] # 2-D
            # Ts_vec[i] = ST.@SVector [orientations[i][1] * T, orientations[i][2] * T, orientations[i][3]] 3-D
        end
    end

    return Ts_vec
end

function state_2_freestream(model::AircraftDynamics, X)
    # calculate freestream
    vx, vy, theta, theta_dot = X[3:6]
    airspeed = [vx, 0.0, -vy] .+ [1e-12, 0.0, 0.0] # in body frame
    vinf = LA.norm(airspeed)
    alpha = asin(airspeed[3] / vinf)
    Omega = ST.@SVector [0.0, theta_dot, 0.0]
    beta = asin(airspeed[2] / vinf)
    freestream = AS.Freestream(vinf, alpha, beta, Omega)

    return freestream
end

function rotor_force(aircraft, parameters, X, freestream, frame; rotors_on = true)
    # extract rotor performance
    Ts = parameters.Ts[:,1]
    orientations = aircraft.rotor_system.orientations
    Ts_vec = body_thrust(Ts, orientations; rotors_on) # in the desired frame, in 3-D

    # calculate forces and moments caused by the rotor in the specified frame
    theta = X[5]
    alpha = freestream.alpha
    Frotor = R_from_body(frame, alpha, theta) * sum(Ts_vec)[[1,3]] * 2 # *2 to mirror
    positions = aircraft.rotor_system.positions
    @assert length(positions) == length(Ts_vec) "Length of rotor positions ($(length(positions))) and thrust vectors ($(length(Ts_vec))) inconsistent."
    # y-axis moment doesn't change between frames
    Mrotor = sum([LA.cross(positions[i] - aircraft.wing_system.system.reference[1].r, Ts_vec[i])[2] for i in 1:length(Ts_vec) ]) # (Ts_vec[2][1] - Ts_vec[1][1]) * (aircraft.wing_system.system.surfaces[2][1].rcp[3] - aircraft.wing_system.system.surfaces[1][1].rcp[3])

    return Frotor, Mrotor
end

"""
    aerodynamics!(model::AircraftDynamics, X, U; rotors_on = true, frame = body::Frame)

Returns x and z forces and y moment in the specified frame.

# Arguments

* `model::AircraftDynamics`
* `X::Vector{Float64}`
* `U::Vector{Float64}`

# Keyword Arguments

* `rotors_on::Bool = true`
* `frame::Frame = Body()`

# Returns

* force
* moment
* power

"""
function aerodynamics!(model::AircraftDynamics, X, U; rotors_on = true, frame::Frame = Body())
    # unpack inputs to AircraftSystems
    aircraft = model.aircraft
    parameters = model.parameters
    environment = model.environment # assumed constant
    stepi = 1 # since there is only 1 step
    stepsymbol = model.stepsymbol

    # unpack states and controls
    x, y, vx, vy, theta, theta_dot, energy = X
    # println("Sherlock!\n\tomegas = $(parameters.omegas)\n\tU = $U")
    parameters.omegas[:,1] .= U # rad/s of the rotors

    # unpack other parameters
    cl_stall = model.cl_stall
    Sref = aircraft.wing_system.system.reference[1].S
    b = aircraft.wing_system.system.reference[1].b
    c = aircraft.wing_system.system.reference[1].c
    g = environment.g

    freestream = state_2_freestream(model, X)

    # solve aerodynamics
    AS.solve_vlm_bem(aircraft, parameters, freestream, environment, (freestream.alpha), stepi, stepsymbol) # solves in the wind axis frame

    # extract aerodynamic forces and moments in wind frame
    CF = ST.@SVector [parameters.CFs[1,1], parameters.CFs[2,1], parameters.CFs[3,1]]
    CM = ST.@SVector [parameters.CMs[1,1], parameters.CMs[2,1], parameters.CMs[3,1]]

    # redimensionalize
    ρ = environment.ρ
    qinf = 0.5 * ρ * freestream.vinf^2

    F = CF .* qinf .* Sref
    M = CM .* (ST.@SVector [b, c, b]) * qinf * Sref

    Qs = parameters.Qs[:,1]
    Ps = Qs .* U
    P = 2 * sum(Ps)

    Frotor, Mrotor = rotor_force(aircraft, parameters, X, freestream, frame; rotors_on)

    # calculate resultants in the specified frame
    force = R_from_wind(frame, freestream.alpha, theta) * F[[1,3]] + Frotor # add thrust
    # force = (ST.@SVector F[[1,3]]) + Frotor # get force vector in x,z plane
    M_y = M[2] + Mrotor

    return force, M_y, P
end

"""
    verbose_aerodynamics!(model::AircraftDynamics, X, U; rotors_on = true, frame::Frame = Body())

Aerodynamics function returns rotor thrust and moment as well.

# Arguments

* `model::AircraftDynamics`
* `X::Vector{Float64}`
* `U::Vector{Float64}`

# Keyword Arguments

* `rotors_on::Bool = true`
* `frame::Frame = Body()`

# Returns

* `force::Vector{Float64}`
* `moment::Float64`
* `power::Float64`
* `rotor_force::Vector{Float64}`
* `rotor_moment::Float64`

"""
# function verbose_aerodynamics!(model::AircraftDynamics, X, U; rotors_on = true, frame::Frame = Body())
#     freestream = state_2_freestream(model, X)
#     total_force, total_moment, power = aerodynamics!(model, X, U; rotors_on, frame)
#     rotor_f, rotor_m = rotor_force(model.aircraft, model.parameters, X, freestream, frame; rotors_on)

#     return total_force, total_moment, power, rotor_f, rotor_m
# end
function verbose_aerodynamics!(model::AircraftDynamics, X, U; rotors_on = true, frame::Frame = Body())
    # unpack inputs to AircraftSystems
    aircraft = model.aircraft
    parameters = model.parameters
    environment = model.environment # assumed constant
    stepi = 1 # since there is only 1 step
    stepsymbol = model.stepsymbol

    # unpack states and controls
    x, y, vx, vy, theta, theta_dot, energy = X
    # println("Sherlock!\n\tomegas = $(parameters.omegas)\n\tU = $U")
    parameters.omegas[:,1] .= U # rad/s of the rotors

    # unpack other parameters
    cl_stall = model.cl_stall
    Sref = aircraft.wing_system.system.reference[1].S
    b = aircraft.wing_system.system.reference[1].b
    c = aircraft.wing_system.system.reference[1].c
    g = environment.g

    freestream = state_2_freestream(model, X)

    # solve aerodynamics
    AS.solve_vlm_bem(aircraft, parameters, freestream, environment, (freestream.alpha), stepi, stepsymbol) # solves in the wind axis frame

    # extract aerodynamic forces and moments in wind frame
    CF = ST.@SVector [parameters.CFs[1,1], parameters.CFs[2,1], parameters.CFs[3,1]]
    CM = ST.@SVector [parameters.CMs[1,1], parameters.CMs[2,1], parameters.CMs[3,1]]

    # redimensionalize
    ρ = environment.ρ
    qinf = 0.5 * ρ * freestream.vinf^2

    F = CF .* qinf .* Sref
    M = CM .* (ST.@SVector [b, c, b]) * qinf * Sref

    Qs = parameters.Qs[:,1]
    Ps = Qs .* U
    P = 2 * sum(Ps)

    Frotor, Mrotor = rotor_force(aircraft, parameters, X, freestream, frame; rotors_on)

    # calculate resultants in the specified frame
    force = R_from_wind(frame, freestream.alpha, theta) * F[[1,3]] + Frotor # add thrust
    # force = (ST.@SVector F[[1,3]]) + Frotor # get force vector in x,z plane
    M_y = M[2] + Mrotor

    println("Sherlock!\n\tparameters.CTs = $(parameters.CTs)\n\tparameters.Ts = $(parameters.Ts)\n\tparameters.omegas = $(parameters.omegas)\n\tvinf = $(freestream.vinf)")
    println("\tC_T = $(parameters.Ts[1] ./ (ρ * (parameters.omegas[1] ./ (2*pi))^2 * (2 * aircraft.rotor_system.rotors[1].Rtip)^4))")

    return force, M_y, P, Frotor, Mrotor
end

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
using Scat
using FFTW
using Test

include("test_filterBank.jl")
include("test_filter.jl")
include("test_scattered.jl")
include("test_scatteringTransform.jl")
include("test_scatter1d.jl")

ext: jl


text:
using KUnet, HDF5, JLD
isdefined(:xtrn) || (@date @load "zn11oparse1.jld")
nc = size(ydev,1)
g0 = 1.0
nbatch = 128
niters = 200
@show net = Layer[KPerceptron(nc, KUnet.kgauss, [g0])]
@date train(net, xdev, ydev; iters=niters, batch=nbatch)
@show (size(net[1].s), nbatch*niters)


text:
module DWDataReader

include("file.jl")

function read(source, kwargs...)
    DWDataReader.File(source; kwargs...)
end

end


text:
struct DeepGPModel{T<:Tuple,L}
    layers::T
    nLayers::Int
    likelihood::L
end

function DeepGPModel(layers...;likelihood=AGP.GaussianLikelihood())
    # Here do some checks with likelihood and last layer
    # Check dimensions between layers
    DeepGPModel{typeof(layers),typeof(likelihood)}(layers,length(layers),likelihood)
end

# Base.:∘(l1::AbstractGPLayer,l2::AbstractGPLayer) = DeepGPModel([l1,l2])
Base.getindex(m::DeepGPModel,i::Int) = m.layers[1]

function ELBO(m::DeepGPModel,X,y)
    expec_log_like(m,X,y) - kl_divergence(m)
end

function kl_divergence(m::DeepGPModel)
    sum(kl_divergence(l) for l in m.layers)
end

propagate(m::DeepGPModel,f) = propagate(Base.tail(m.layers),propagate(first(m.layers),f))
propagate(::Tuple{}, f) = f
propagate(layers::Tuple,f) = propagate(Base.tail(layers),propagate(first(layers),f))

function expec_log_like(m,X,y,nSamples=1)
    f = [X for _ in 1:nSamples]
    f_final = propagate(m,f)
    # for l in m.layers
    #     f = propagate(l,f)
    # end
    loss = 0.0
    for i in 1:nSamples
        loss += logpdf(m,f_final[i],y)/nSamples
    end
    return loss
end

function logpdf(m::DeepGPModel,f,y)
    tot = 0.0
    for i in 1:length(y)
        tot += AGP.logpdf(m.likelihood,f[i],y[i]) #Would not work for multiclass for example
    end
    return tot
end


text:
@printf("\nTesting linear operators calculus rules\n")

##########################
##### test Compose #######
##########################
m1, m2, m3 = 4, 7, 3
A1 = randn(m2, m1)
A2 = randn(m3, m2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)

opC = Compose(opA2,opA1)
x = randn(m1)
y1 = test_op(opC, x, randn(m3), verb)
y2 = A2*A1*x
@test all(norm.(y1 .- y2) .<= 1e-12)

# test Compose longer
m1, m2, m3, m4 = 4, 7, 3, 2
A1 = randn(m2, m1)
A2 = randn(m3, m2)
A3 = randn(m4, m3)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)

opC1 = Compose(opA3,Compose(opA2,opA1))
opC2 = Compose(Compose(opA3,opA2),opA1)
x = randn(m1)
y1 = test_op(opC1, x, randn(m4), verb)
y2 = test_op(opC2, x, randn(m4), verb)
y3 = A3*A2*A1*x
@test all(norm.(y1 .- y2) .<= 1e-12)
@test all(norm.(y3 .- y2) .<= 1e-12)

#test Compose special cases
@test typeof(opA1*Eye(m1)) == typeof(opA1)
@test typeof(Eye(m2)*opA1) == typeof(opA1)
@test typeof(Eye(m2)*Eye(m2)) == typeof(Eye(m2))

opS1 = Compose(opA2,opA1)
opS1c = Scale(pi,opS1)
@test typeof(opS1c.A[end]) <: Scale

#properties
@test is_sliced(opC)            == false
@test is_linear(opC1)           == true
@test is_null(opC1)             == false
@test is_eye(opC1)              == false
@test is_diagonal(opC1)         == false
@test is_AcA_diagonal(opC1)     == false
@test is_AAc_diagonal(opC1)     == false
@test is_orthogonal(opC1)       == false
@test is_invertible(opC1)       == false
@test is_full_row_rank(opC1)    == false
@test is_full_column_rank(opC1) == false

# properties special case
opC = DCT((5,))*GetIndex((10,), 1:5)
@test is_sliced(opC)           == true
@test is_AAc_diagonal(opC)     == true
@test diag_AAc(opC)            == 1.0

d = randn(5)
opC = DiagOp(d)*GetIndex((10,), 1:5)
@test is_sliced(opC)           == true
@test is_diagonal(opC)         == true
@test diag(opC)                == d

# displacement test
m1, m2, m3, m4, m5 = 4, 7, 3, 2, 11
A1 = randn(m2, m1)
A2 = randn(m3, m2)
A3 = randn(m4, m3)
A4 = randn(m5, m4)
d1 = randn(m2)
d2 = pi
d3 = 0.0
d4 = randn(m5)
opA1 = AffineAdd(MatrixOp(A1),d1)
opA2 = AffineAdd(MatrixOp(A2),d2)
opA3 = MatrixOp(A3)
opA4 = AffineAdd(MatrixOp(A4),d4,false)

opC  = Compose(Compose(Compose(opA4,opA3),opA2),opA1)

x = randn(m1)

@test norm( opC*x - (A4*(A3*( A2*( A1*x+d1 ) .+ d2 ) .+ d3)-d4) )  < 1e-9
@test norm( displacement(opC) - ( A4*(A3*(A2*d1 .+d2 ) .+ d3)-d4) ) < 1e-9

opA4 = MatrixOp(A4)
opC  = AffineAdd(Compose(Compose(Compose(opA4,opA3),opA2),opA1),d4)
@test norm( opC*x - (A4*(A3*( A2*( A1*x+d1 ) .+ d2 ) .+ d3)+d4) )  < 1e-9
@test norm( displacement(opC) - ( A4*(A3*(A2*d1 .+ d2) .+ d3)+d4) ) < 1e-9

@test norm( remove_displacement(opC)*x - (A4*(A3*( A2*( A1*x ) ) ) ) )  < 1e-9

########################
### test DCAT    #######
########################

m1, n1, m2, n2 = 4, 7, 5, 2
A1 = randn(m1, n1)
A2 = randn(m2, n2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opD = DCAT(opA1, opA2)
x1 = randn(n1)
x2 = randn(n2)
y1 = test_op(opD, ArrayPartition(x1, x2), ArrayPartition(randn(m1),randn(m2)), verb)
y2 = ArrayPartition(A1*x1, A2*x2)
@test norm(y1 .- y2) .<= 1e-12

# test DCAT longer

m1, n1, m2, n2, m3, n3 = 4, 7, 5, 2, 5, 5
A1 = randn(m1, n1)
A2 = randn(m2, n2)
A3 = randn(m3, n3)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
opD = DCAT(opA1, opA2, opA3)
x1 = randn(n1)
x2 = randn(n2)
x3 = randn(n3)
y1 = test_op(opD, ArrayPartition(x1, x2, x3), ArrayPartition(randn(m1),randn(m2),randn(m3)), verb)
y2 = ArrayPartition(A1*x1, A2*x2, A3*x3)
@test norm(y1 .- y2) .<= 1e-12

#properties
@test is_linear(opD)           == true
@test is_null(opD)             == false
@test is_eye(opD)              == false
@test is_diagonal(opD)         == false
@test is_AcA_diagonal(opD)     == false
@test is_AAc_diagonal(opD)     == false
@test is_orthogonal(opD)       == false
@test is_invertible(opD)       == false
@test is_full_row_rank(opD)    == false
@test is_full_column_rank(opD) == false

# DCAT of Eye

n1, n2 = 4, 7
x1 = randn(n1)
x2 = randn(n2)

opD = Eye(ArrayPartition(x1,x2))
y1 = test_op(opD, ArrayPartition(x1, x2), ArrayPartition(randn(n1),randn(n2)), verb)

#properties
@test is_linear(opD)           == true
@test is_null(opD)             == false
@test is_eye(opD)              == true
@test is_diagonal(opD)         == true
@test is_AcA_diagonal(opD)     == true
@test is_AAc_diagonal(opD)     == true
@test is_orthogonal(opD)       == true
@test is_invertible(opD)       == true
@test is_full_row_rank(opD)    == true
@test is_full_column_rank(opD) == true

@test diag(opD) == 1
@test diag_AcA(opD) == 1
@test diag_AAc(opD) == 1

# displacement DCAT

m1, n1, m2, n2, m3, n3 = 4, 7, 5, 2, 5, 5
A1 = randn(m1, n1)
A2 = randn(m2, n2)
A3 = randn(m3, n3)
d1 = randn(m1)
d2 = randn(m2)
d3 = randn(m3)
opA1 = AffineAdd(MatrixOp(A1), d1)
opA2 = AffineAdd(MatrixOp(A2), d2)
opA3 = AffineAdd(MatrixOp(A3), d3)
opD = DCAT(opA1, opA2, opA3)
x1 = randn(n1)
x2 = randn(n2)
x3 = randn(n3)
y1 = opD*ArrayPartition(x1,x2,x3)
y2 = ArrayPartition(A1*x1+d1, A2*x2+d2, A3*x3+d3)
@test norm(y1 .- y2) .<= 1e-12
@test norm(displacement(opD) .- ArrayPartition(d1,d2,d3)) .<= 1e-12
y1 = remove_displacement(opD)*ArrayPartition(x1,x2,x3)
y2 = ArrayPartition(A1*x1, A2*x2, A3*x3)
@test norm(y1 .- y2) .<= 1e-12

#######################
## test HCAT    #######
#######################

m, n1, n2 = 4, 7, 5
A1 = randn(m, n1)
A2 = randn(m, n2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opH = HCAT(opA1, opA2)
x1 = randn(n1)
x2 = randn(n2)
y1 = test_op(opH, ArrayPartition(x1, x2), randn(m), verb)
y2 = A1*x1 + A2*x2
@test norm(y1-y2) <= 1e-12

#permuatation
p = [2;1]
opHp = opH[p]
y1 = test_op(opHp, ArrayPartition(x2, x1), randn(m), verb)
@test norm(y1-y2) <= 1e-12

# test HCAT longer

m, n1, n2, n3 = 4, 7, 5, 6
A1 = randn(m, n1)
A2 = randn(m, n2)
A3 = randn(m, n3)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
opH = HCAT(opA1, opA2, opA3)
x1 = randn(n1)
x2 = randn(n2)
x3 = randn(n3)
y1 = test_op(opH, ArrayPartition(x1, x2, x3), randn(m), verb)
y2 = A1*x1 + A2*x2 + A3*x3
@test norm(y1-y2) <= 1e-12

# test HCAT of HCAT
opHH = HCAT(opH, opA2, opA3)
y1 = test_op(opHH, ArrayPartition(x1, x2, x3, x2, x3), randn(m), verb)
y2 = A1*x1 + A2*x2 + A3*x3 + A2*x2 + A3*x3
@test norm(y1-y2) <= 1e-12

opHH = HCAT(opH, opH, opA3)
x = ArrayPartition(x1, x2, x3, x1, x2, x3, x3)
y1 = test_op(opHH, x, randn(m), verb)
y2 = A1*x1 + A2*x2 + A3*x3 + A1*x1 + A2*x2 + A3*x3 + A3*x3
@test norm(y1-y2) <= 1e-12

opA3 = MatrixOp(randn(n1,n1))
@test_throws Exception HCAT(opA1,opA2,opA3)
opF = AbstractOperators.DFT(Complex{Float64},(m,))
@test_throws Exception HCAT(opA1,opF,opA2)

# test utilities

# permutation
p = randperm(ndoms(opHH,2))
opHP = AbstractOperators.permute(opHH,p)

xp = ArrayPartition(x.x[p]...)

y1 = test_op(opHP, xp, randn(m), verb)

pp = randperm(ndoms(opHH,2))
opHPP = AbstractOperators.permute(opHH,pp)
xpp = ArrayPartition(x.x[pp]...)
y1 = test_op(opHPP, xpp, randn(m), verb)

#properties
m, n1, n2, n3 = 4, 7, 5, 6
A1 = randn(m, n1)
A2 = randn(m, n2)
A3 = randn(m, n3)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
op = HCAT(opA1, opA2, opA3)
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == false

d = randn(n1).+im.*randn(n1)
op = HCAT(DiagOp(d), AbstractOperators.DFT(Complex{Float64},n1))
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == false

@test diag_AAc(op) == d .* conj(d) .+ n1

y1 = randn(n1) .+ im .* randn(n1)
@test norm(op*(op'*y1).-diag_AAc(op).*y1) <1e-12

#test displacement

m, n1, n2 = 4, 7, 5
A1 = randn(m, n1)
A2 = randn(m, n2)
d1 = randn(m)
d2 = randn(m)
opA1 = AffineAdd(MatrixOp(A1), d1)
opA2 = AffineAdd(MatrixOp(A2), d2)
opH = HCAT(opA1, opA2)
x1 = randn(n1)
x2 = randn(n2)
y1 = opH*ArrayPartition(x1,x2)
y2 = A1*x1+d1 + A2*x2+d2
@test norm(y1-y2) <= 1e-12
y1 = remove_displacement(opH)*ArrayPartition(x1,x2)
y2 = A1*x1 + A2*x2
@test norm(y1-y2) <= 1e-12

########################
### test Reshape #######
########################

m, n = 8, 4
dim_out = (2, 2, 2)
A1 = randn(m, n)
opA1 = MatrixOp(A1)
opR = Reshape(opA1, dim_out)
opR = Reshape(opA1, dim_out...)
x1 = randn(n)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = reshape(A1*x1, dim_out)
@test norm(y1-y2) <= 1e-12

@test_throws Exception Reshape(opA1,(2,2,1))

@test is_null(opR)             == is_null(opA1)
@test is_eye(opR)              == is_eye(opA1)
@test is_diagonal(opR)         == is_diagonal(opA1)
@test is_AcA_diagonal(opR)     == is_AcA_diagonal(opA1)
@test is_AAc_diagonal(opR)     == is_AAc_diagonal(opA1)
@test is_orthogonal(opR)       == is_orthogonal(opA1)
@test is_invertible(opR)       == is_invertible(opA1)
@test is_full_row_rank(opR)    == is_full_row_rank(opA1)
@test is_full_column_rank(opR) == is_full_column_rank(opA1)

# testing displacement
m, n = 8, 4
dim_out = (2, 2, 2)
A1 = randn(m, n)
d1 = randn(m)
opA1 = AffineAdd(MatrixOp(A1),d1)
opR = Reshape(opA1, dim_out)
x1 = randn(n)
y1 = opR*x1
y2 = reshape(A1*x1+d1, dim_out)
@test norm(y1-y2) <= 1e-12
y1 = remove_displacement(opR)*x1
y2 = reshape(A1*x1, dim_out)
@test norm(y1-y2) <= 1e-12

#######################
## test Scale   #######
#######################

m, n = 8, 4
coeff = pi
A1 = randn(m, n)
opA1 = MatrixOp(A1)
opS = Scale(coeff, opA1)
x1 = randn(n)
y1 = test_op(opS, x1, randn(m), verb)
y2 = coeff*A1*x1
@test norm(y1-y2) <= 1e-12

coeff2 = 3
opS2 = Scale(coeff2, opS)
y1 = test_op(opS2, x1, randn(m), verb)
y2 = coeff2*coeff*A1*x1
@test norm(y1-y2) <= 1e-12

opF = AbstractOperators.DFT(m,n)
opS = Scale(coeff, opF)
x1 = randn(m,n)
y1 = test_op(opS, x1, fft(randn(m,n)), verb)
y2 = coeff*(fft(x1))
@test norm(y1-y2) <= 1e-12

opS = Scale(coeff, opA1)
@test is_null(opS)             == is_null(opA1)
@test is_eye(opS)              == is_eye(opA1)
@test is_diagonal(opS)         == is_diagonal(opA1)
@test is_AcA_diagonal(opS)     == is_AcA_diagonal(opA1)
@test is_AAc_diagonal(opS)     == is_AAc_diagonal(opA1)
@test is_orthogonal(opS)       == is_orthogonal(opA1)
@test is_invertible(opS)       == is_invertible(opA1)
@test is_full_row_rank(opS)    == is_full_row_rank(opA1)
@test is_full_column_rank(opS) == is_full_column_rank(opA1)

op = Scale(-4.0,AbstractOperators.DFT(10))
@test is_AAc_diagonal(op)     == true
@test diag_AAc(op) == 16*10

op = Scale(-4.0,ZeroPad((10,), 20))
@test is_AcA_diagonal(op)     == true
@test diag_AcA(op) == 16

# special case, Scale of DiagOp gets a DiagOp
d = randn(10)
op = Scale(3,DiagOp(d))
@test typeof(op) <: DiagOp
@test norm(diag(op) - 3 .*d) < 1e-12

# Scale with imaginary coeff gives error
m, n = 8, 4
coeff = im
A1 = randn(m, n)
opA1 = MatrixOp(A1)
@test_throws ErrorException Scale(coeff, opA1)

## testing displacement
m, n = 8, 4
coeff = pi
A1 = randn(m, n)
d1 = randn(m)
opA1 = AffineAdd(MatrixOp(A1),d1)
opS = Scale(coeff, opA1)
x1 = randn(n)
y1 = opS*x1
y2 = coeff*(A1*x1+d1)
@test norm(y1-y2) <= 1e-12
y1 = remove_displacement(opS)*x1
y2 = coeff*(A1*x1)
@test norm(y1-y2) <= 1e-12

########################
### test Sum     #######
########################

m,n = 5,7
A1 = randn(m,n)
A2 = randn(m,n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opS = Sum(opA1,opA2)
x1 = randn(n)
y1 = test_op(opS, x1, randn(m), verb)
y2 = A1*x1+A2*x1
@test norm(y1-y2) <= 1e-12

#test Sum longer
m,n = 5,7
A1 = randn(m,n)
A2 = randn(m,n)
A3 = randn(m,n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
opS = Sum(opA1,opA2,opA3)
x1 = randn(n)
y1 = test_op(opS, x1, randn(m), verb)
y2 = A1*x1+A2*x1+A3*x1
@test norm(y1-y2) <= 1e-12

opA3 = MatrixOp(randn(m,m))
@test_throws Exception Sum(opA1,opA3)
opF = AbstractOperators.DFT(Float64,(m,))
@test_throws Exception Sum(opF,opA3)

@test is_null(opS)             == false
@test is_eye(opS)              == false
@test is_diagonal(opS)         == false
@test is_AcA_diagonal(opS)     == false
@test is_AAc_diagonal(opS)     == false
@test is_orthogonal(opS)       == false
@test is_invertible(opS)       == false
@test is_full_row_rank(opS)    == true
@test is_full_column_rank(opS) == false

d = randn(10)
op = Sum(Scale(-3.1,Eye(10)),DiagOp(d))
@test is_diagonal(op)         == true
@test norm(   diag(op) - (d .-3.1)  )<1e-12

#test displacement of sum
m,n = 5,7
A1 = randn(m,n)
A2 = randn(m,n)
A3 = randn(m,n)
d1 = randn(m)
d2 = pi
d3 = randn(m)
opA1 = AffineAdd(MatrixOp(A1), d1)
opA2 = AffineAdd(MatrixOp(A2), d2)
opA3 = AffineAdd(MatrixOp(A3), d3)
opS = Sum(opA1,opA2,opA3)
x1 = randn(n)
y2 = A1*x1+A2*x1+A3*x1+d1.+d2+d3
@test norm(opS*x1-y2) <= 1e-12
@test norm(displacement(opS) - (d1.+d2+d3)) <= 1e-12
y2 = A1*x1+A2*x1+A3*x1
@test norm(remove_displacement(opS)*x1-y2) <= 1e-12

###################################
###### test AdjointOperator #######
###################################

m,n = 5,7
A1 = randn(m,n)
opA1 = MatrixOp(A1)
opA1t = MatrixOp(A1')
opT = AdjointOperator(opA1)
x1 = randn(m)
y1 = test_op(opT, x1, randn(n), verb)
y2 = A1'*x1
@test norm(y1-y2) <= 1e-12

@test is_null(opT)             == is_null(opA1t)
@test is_eye(opT)              == is_eye(opA1t)
@test is_diagonal(opT)         == is_diagonal(opA1t)
@test is_AcA_diagonal(opT)     == is_AcA_diagonal(opA1t)
@test is_AAc_diagonal(opT)     == is_AAc_diagonal(opA1t)
@test is_orthogonal(opT)       == is_orthogonal(opA1t)
@test is_invertible(opT)       == is_invertible(opA1t)
@test is_full_row_rank(opT)    == is_full_row_rank(opA1t)
@test is_full_column_rank(opT) == is_full_column_rank(opA1t)

d = randn(3)
op = AdjointOperator(DiagOp(d))
@test is_diagonal(op) == true
@test diag(op) == d

op = AdjointOperator(ZeroPad((10,),5))
@test is_AcA_diagonal(op) == false
@test is_AAc_diagonal(op) == true
@test diag_AAc(op) == 1

#############################
######## test VCAT    #######
#############################

m1, m2, n = 4, 7, 5
A1 = randn(m1, n)
A2 = randn(m2, n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opV = VCAT(opA1, opA2)
x1 = randn(n)
y1 = test_op(opV, x1, ArrayPartition(randn(m1), randn(m2)), verb)
y2 = ArrayPartition(A1*x1, A2*x1)
@test norm(y1 - y2) .<= 1e-12

m1, n = 4, 5
A1 = randn(m1, n)+im*randn(m1, n)
opA1 = MatrixOp(A1)
opA2 = AbstractOperators.DFT(n)'
opV = VCAT(opA1, opA2)
x1 = randn(n)+im*randn(n)
y1 = test_op(opV, x1, ArrayPartition(randn(m1)+im*randn(m1), randn(n)), verb)
y2 = ArrayPartition(A1*x1, opA2*x1)
@test norm(y1 - y2) .<= 1e-12

#test VCAT longer
m1, m2, m3, n = 4, 7, 3, 5
A1 = randn(m1, n)
A2 = randn(m2, n)
A3 = randn(m3, n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
opV = VCAT(opA1, opA2, opA3)
x1 = randn(n)
y1 = test_op(opV, x1, ArrayPartition(randn(m1), randn(m2), randn(m3)), verb)
y2 = ArrayPartition(A1*x1, A2*x1, A3*x1)
@test norm(y1 - y2) .<= 1e-12

#test VCAT of VCAT
opVV = VCAT(opV,opA3)
y1 = test_op(opVV, x1, ArrayPartition(randn(m1), randn(m2), randn(m3), randn(m3)), verb)
y2 = ArrayPartition(A1*x1, A2*x1, A3*x1, A3*x1)
@test norm(y1 .- y2) <= 1e-12

opVV = VCAT(opA1,opV,opA3)
y1 = test_op(opVV, x1, ArrayPartition(randn(m1), randn(m1), randn(m2), randn(m3), randn(m3)), verb)
y2 = ArrayPartition(A1*x1, A1*x1, A2*x1, A3*x1, A3*x1)
@test norm(y1 - y2) <= 1e-12

opA3 = MatrixOp(randn(m1,m1))
@test_throws Exception VCAT(opA1,opA2,opA3)
opF = AbstractOperators.DFT(Complex{Float64},(n,))
@test_throws Exception VCAT(opA1,opF,opA2)

###properties
m1, m2, m3, n = 4, 7, 3, 5
A1 = randn(m1, n)
A2 = randn(m2, n)
A3 = randn(m3, n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
op = VCAT(opA1, opA2, opA3)
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == false
@test is_full_column_rank(op) == true

op = VCAT(AbstractOperators.DFT(Complex{Float64},10), Eye(Complex{Float64},10) )
@test is_AcA_diagonal(op)     == true
@test diag_AcA(op) == 11

##test displacement
m1, m2, n = 4, 7, 5
A1 = randn(m1, n)
A2 = randn(m2, n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
d1 = randn(m1)
d2 = randn(m2)
opV = VCAT(AffineAdd(opA1,d1), AffineAdd(opA2,d2))
x1 = randn(n)
y1 = opV*x1
y2 = ArrayPartition(A1*x1+d1, A2*x1+d2)
@test norm(y1 - y2) <= 1e-12
y1 = remove_displacement(opV)*x1
y2 = ArrayPartition(A1*x1, A2*x1)
@test norm(y1 - y2) <= 1e-12

#########################
#### test BroadCast #####
#########################

m, n = 8, 4
dim_out = (m, 10)
A1 = randn(m, n)
opA1 = MatrixOp(A1)
opR = BroadCast(opA1, dim_out)
x1 = randn(n)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= A1*x1
@test norm(y1-y2) <= 1e-12

m, n, l, k = 8, 4, 5, 7
dim_out = (m, n, l, k)
opA1 = Eye(m,n)
opR = BroadCast(opA1, dim_out)
x1 = randn(m,n)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= x1
@test norm(y1-y2) <= 1e-12

@test_throws Exception BroadCast(opA1,(m,m))

m, l = 1, 5
dim_out = (m, l)
opA1 = Eye(m)
opR = BroadCast(opA1, dim_out)
x1 = randn(m)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= x1
@test norm(y1-y2) <= 1e-12

#colum in - matrix out
m, l = 4, 5
dim_out = (m, l)
opA1 = Eye(1,l)
opR = BroadCast(opA1, dim_out)
x1 = randn(1,l)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= x1
@test norm(y1-y2) <= 1e-12

op = HCAT(Eye(m,l),opR)
x1 = ArrayPartition(randn(m,l),randn(1,l))
y1 = test_op(op, x1, randn(dim_out), verb)
y2 = x1.x[1].+x1.x[2]
@test norm(y1-y2) <= 1e-12

m, n, l  = 2, 5, 8
dim_out = (m, n, l)
opA1 = Eye(m)
opR = BroadCast(opA1, dim_out)
x1 = randn(m)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= x1
@test norm(y1-y2) <= 1e-12

m, n, l  = 1, 5, 8
dim_out = (m, n, l)
opA1 = Eye(m)
opR = BroadCast(opA1, dim_out)
x1 = randn(m)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= x1
@test norm(y1-y2) <= 1e-12

m, n, l  = 1, 5, 8
dim_out = (m, n, l)
opA1 = Scale(2.4,Eye(m))
opR = BroadCast(opA1, dim_out)
x1 = randn(m)
y1 = test_op(opR, x1, randn(dim_out), verb)
y2 = zeros(dim_out)
y2 .= 2.4*x1
@test norm(y1-y2) <= 1e-12

@test is_null(opR)             == is_null(opA1)
@test is_eye(opR)              == false
@test is_diagonal(opR)         == false
@test is_AcA_diagonal(opR)     == false
@test is_AAc_diagonal(opR)     == false
@test is_orthogonal(opR)       == false
@test is_invertible(opR)       == false
@test is_full_row_rank(opR)    == false
@test is_full_column_rank(opR) == false

# test displacement

m, n = 8, 4
dim_out = (m, 10)
A1 = randn(m, n)
d1 = randn(m)
opA1 = AffineAdd(MatrixOp(A1), d1)
opR = BroadCast(opA1, dim_out)
x1 = randn(n)
y1 = opR*x1
y2 = zeros(dim_out)
y2 .= A1*x1+d1
@test norm(y1-y2) <= 1e-12
x1 = randn(n)
y1 = remove_displacement(opR)*x1
y2 = zeros(dim_out)
y2 .= A1*x1
@test norm(y1-y2) <= 1e-12

###########################
##### test AffineAdd  #####
###########################

n,m = 5,6
A = randn(n,m)
opA = MatrixOp(A)
d = randn(n)
T = AffineAdd(opA,d)

println(T)
x1 = randn(m)
y1 = T*x1
@test norm(y1-(A*x1+d)) <1e-9
r = randn(n)
@test norm(T'*r-(A'*r)) <1e-9
@test displacement(T) == d
@test norm(remove_displacement(T)*x1-A*x1) <1e-9

# with sign
T = AffineAdd(opA,d,false)
@test sign(T) == -1

println(T)
x1 = randn(m)
y1 = T*x1
@test norm(y1-(A*x1-d)) <1e-9
r = randn(n)
@test norm(T'*r-(A'*r)) <1e-9
@test displacement(T) == -d
@test norm(remove_displacement(T)*x1-A*x1) <1e-9

# with scalar
T = AffineAdd(opA,pi)
@test sign(T) == 1

println(T)
x1 = randn(m)
y1 = T*x1
@test norm(y1-(A*x1.+pi)) <1e-9
r = randn(n)
@test norm(T'*r-(A'*r)) < 1e-9
@test displacement(T) .- pi < 1e-9
@test norm(remove_displacement(T)*x1-A*x1) <1e-9

@test_throws DimensionMismatch AffineAdd(MatrixOp(randn(2,5)),randn(5))
@test_throws ErrorException AffineAdd(AbstractOperators.DFT(4),randn(4))
AffineAdd(AbstractOperators.DFT(4),pi)
@test_throws ErrorException AffineAdd(Eye(4),im*pi)

# with scalar and vector
d = randn(n)
T = AffineAdd(AffineAdd(opA,pi),d,false)

println(T)
x1 = randn(m)
y1 = T*x1
@test norm(y1-(A*x1 .+ pi .- d)) <1e-9
r = randn(n)
@test norm(T'*r-(A'*r)) < 1e-9
@test norm(displacement(T) .- (pi .-d )) < 1e-9

T2 = remove_displacement(T)
@test norm(T2*x1-(A*x1)) <1e-9

# permute AddAffine
n,m = 5,6
A = randn(n,m)
d = randn(n)
opH = HCAT(Eye(n),MatrixOp(A))
x = ArrayPartition(randn(n),randn(m))
opHT = AffineAdd(opH,d)

@test norm(opHT*x-(x.x[1]+A*x.x[2].+d)) < 1e-12
p = [2;1]
@test norm(AbstractOperators.permute(opHT,p)*ArrayPartition(x.x[p]...)-(x.x[1]+A*x.x[2].+d)) < 1e-12

############################
####### test combin. #######
############################

## test Compose of HCAT
m1, m2, m3, m4 = 4, 7, 3, 2
A1 = randn(m3, m1)
A2 = randn(m3, m2)
A3 = randn(m4, m3)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
opH = HCAT(opA1,opA2)
opC = Compose(opA3,opH)
x1, x2 = randn(m1), randn(m2)
y1 = test_op(opC, ArrayPartition(x1,x2), randn(m4), verb)

y2 = A3*(A1*x1+A2*x2)

@test norm(y1-y2) < 1e-9

opCp = AbstractOperators.permute(opC,[2,1])
y1 = test_op(opCp, ArrayPartition(x2,x1), randn(m4), verb)

@test norm(y1-y2) < 1e-9

## test HCAT of Compose of HCAT
m5 = 10
A4 = randn(m4,m5)
x3 = randn(m5)
opHC = HCAT(opC,MatrixOp(A4))
x = ArrayPartition(x1,x2,x3)
y1 = test_op(opHC, x, randn(m4), verb)

@test norm(y1-(y2+A4*x3)) < 1e-9

p = randperm(ndoms(opHC,2))
opHP = AbstractOperators.permute(opHC,p)

xp = ArrayPartition(x.x[p]...)

y1 = test_op(opHP, xp, randn(m4), verb)

pp = randperm(ndoms(opHC,2))
opHPP = AbstractOperators.permute(opHC,pp)
xpp = ArrayPartition(x.x[pp]...)
y1 = test_op(opHPP, xpp, randn(m4), verb)

# test VCAT of HCAT's
m1, m2, n1 = 4, 7, 3
A1 = randn(n1, m1)
A2 = randn(n1, m2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opH1 = HCAT(opA1,opA2)

m1, m2, n2 = 4, 7, 5
A3 = randn(n2, m1)
A4 = randn(n2, m2)
opA3 = MatrixOp(A3)
opA4 = MatrixOp(A4)
opH2 = HCAT(opA3,opA4)

opV = VCAT(opH1,opH2)
x1, x2 = randn(m1), randn(m2)
y1 = test_op(opV, ArrayPartition(x1,x2), ArrayPartition(randn(n1),randn(n2)), verb)
y2 = ArrayPartition(A1*x1+A2*x2,A3*x1+A4*x2)
@test norm(y1 - y2) <= 1e-12

# test VCAT of HCAT's with complex num
m1, m2, n1 = 4, 7, 5
A1 = randn(n1, m1)+im*randn(n1, m1)
opA1 = MatrixOp(A1)
opA2 = AbstractOperators.DFT(n1)
opH1 = HCAT(opA1,opA2)

m1, m2, n2 = 4, 7, 5
A3 = randn(n2, m1)+im*randn(n2,m1)
opA3 = MatrixOp(A3)
opA4 = AbstractOperators.DFT(n2)
opH2 = HCAT(opA3,opA4)

opV = VCAT(opH1,opH2)
x1, x2 = randn(m1)+im*randn(m1), randn(n2)
y1 = test_op(opV, ArrayPartition(x1,x2), ArrayPartition(randn(n1)+im*randn(n1),randn(n2)+im*randn(n2)), verb)
y2 = ArrayPartition(A1*x1+fft(x2),A3*x1+fft(x2))
@test norm(y1 - y2) <= 1e-12

# test HCAT of VCAT's

n1, n2, m1, m2 = 3, 5, 4, 7
A = randn(m1, n1); opA = MatrixOp(A)
B = randn(m1, n2); opB = MatrixOp(B)
C = randn(m2, n1); opC = MatrixOp(C)
D = randn(m2, n2); opD = MatrixOp(D)
opV = HCAT(VCAT(opA, opC), VCAT(opB, opD))
x1 = randn(n1)
x2 = randn(n2)
y1 = test_op(opV, ArrayPartition(x1, x2), ArrayPartition(randn(m1), randn(m2)), verb)
y2 = ArrayPartition(A*x1 + B*x2, C*x1 + D*x2)

@test norm(y1 - y2) <= 1e-12

# test Sum of HCAT's

m, n1, n2, n3 = 4, 7, 5, 3
A1 = randn(m, n1)
A2 = randn(m, n2)
A3 = randn(m, n3)
B1 = randn(m, n1)
B2 = randn(m, n2)
B3 = randn(m, n3)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA3 = MatrixOp(A3)
opB1 = MatrixOp(B1)
opB2 = MatrixOp(B2)
opB3 = MatrixOp(B3)
opHA = HCAT(opA1, opA2, opA3)
opHB = HCAT(opB1, opB2, opB3)
opS = Sum(opHA, opHB)
x1 = randn(n1)
x2 = randn(n2)
x3 = randn(n3)
y1 = test_op(opS, ArrayPartition(x1, x2, x3), randn(m), verb)
y2 = A1*x1 + B1*x1 + A2*x2 + B2*x2 + A3*x3 + B3*x3

@test norm(y1-y2) <= 1e-12

p = [3;2;1]
opSp = AbstractOperators.permute(opS,p)
y1 = test_op(opSp, ArrayPartition(((x1, x2, x3)[p])...), randn(m), verb)

# test Sum of VCAT's

m1, m2, n = 4, 7, 5
A1 = randn(m1, n)
A2 = randn(m2, n)
B1 = randn(m1, n)
B2 = randn(m2, n)
C1 = randn(m1, n)
C2 = randn(m2, n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opB1 = MatrixOp(B1)
opB2 = MatrixOp(B2)
opC1 = MatrixOp(C1)
opC2 = MatrixOp(C2)
opVA = VCAT(opA1, opA2)
opVB = VCAT(opB1, opB2)
opVC = VCAT(opC1, opC2)
opS = Sum(opVA, opVB, opVC)
x = randn(n)
y1 = test_op(opS, x, ArrayPartition(randn(m1), randn(m2)), verb)
y2 = ArrayPartition(A1*x + B1*x +C1*x, A2*x + B2*x + C2*x)

@test norm(y1 - y2) .<= 1e-12

# test Scale of DCAT

m1, n1 = 4, 7
m2, n2 = 3, 5
A1 = randn(m1, n1)
A2 = randn(m2, n2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opD = DCAT(opA1, opA2)
coeff = randn()
opS = Scale(coeff, opD)
x1 = randn(n1)
x2 = randn(n2)
y = test_op(opS, ArrayPartition(x1, x2), ArrayPartition(randn(m1), randn(m2)), verb)
z = ArrayPartition(coeff*A1*x1, coeff*A2*x2)

@test norm(y - z) <= 1e-12

# test Scale of VCAT

m1, m2, n = 4, 3, 7
A1 = randn(m1, n)
A2 = randn(m2, n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opV = VCAT(opA1, opA2)
coeff = randn()
opS = Scale(coeff, opV)
x = randn(n)
y = test_op(opS, x, ArrayPartition(randn(m1), randn(m2)), verb)
z = ArrayPartition(coeff*A1*x, coeff*A2*x)

@test norm(y - z) <= 1e-12

# test Scale of HCAT

m, n1, n2 = 4, 3, 7
A1 = randn(m, n1)
A2 = randn(m, n2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opH = HCAT(opA1, opA2)
coeff = randn()
opS = Scale(coeff, opH)
x1 = randn(n1)
x2 = randn(n2)
y = test_op(opS, ArrayPartition(x1, x2), randn(m), verb)
z = coeff*(A1*x1 + A2*x2)

@test norm(y - z) <= 1e-12

# test DCAT of HCATs

m1, m2, n1, n2, l1, l2, l3 = 2,3,4,5,6,7,8
A1 = randn(m1, n1)
A2 = randn(m1, n2)
B1 = randn(m2, n1)
B2 = randn(m2, n2)
B3 = randn(m2, n2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opH1 = HCAT(opA1, opA2)
opB1 = MatrixOp(B1)
opB2 = MatrixOp(B2)
opB3 = MatrixOp(B3)
opH2 = HCAT(opB1, opB2, opB3)

op = DCAT(opA1, opH2)
x =  ArrayPartition(randn.(size(op,2))...)
y0 = ArrayPartition(randn.(size(op,1))...)
y = test_op(op, x, y0, verb)

op = DCAT(opH1, opH2)
x =  ArrayPartition(randn.(size(op,2))...)
y0 = ArrayPartition(randn.(size(op,1))...)
y = test_op(op, x, y0, verb)

p = randperm(ndoms(op,2))
y2 = op[p]*ArrayPartition(x.x[p]...)

@test norm(y - y2) <= 1e-8

# test Scale of Sum

m,n = 5,7
A1 = randn(m,n)
A2 = randn(m,n)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opS = Sum(opA1,opA2)
coeff = pi
opSS = Scale(coeff,opS)
x1 = randn(n)
y1 = test_op(opSS, x1, randn(m), verb)
y2 = coeff*(A1*x1+A2*x1)
@test norm(y1-y2) <= 1e-12

# test Scale of Compose

m1, m2, m3 = 4, 7, 3
A1 = randn(m2, m1)
A2 = randn(m3, m2)
opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)

coeff = pi
opC = Compose(opA2,opA1)
opS = Scale(coeff,opC)
x = randn(m1)
y1 = test_op(opS, x, randn(m3), verb)
y2 = coeff*(A2*A1*x)
@test all(norm.(y1 .- y2) .<= 1e-12)


text:
using Test

for file in readdir(@__DIR__)
    if file in ["Utilities.jl"]
        continue
    end
    @testset "$(file)" begin
        include(file)
    end
end


text:
module IL_simulations

  using Parameters
  using DelimitedFiles
  using PDBTools

  include("./mdp_files.jl")
  include("dict_data.jl")
  include("./system_data.jl")
  include("./conc_calc.jl")
  include("./measure_protein.jl")
  include("./write_pack_input.jl")
  include("./topology.jl")
  include("./mddf_calculations.jl")
  include("./util.jl")

end


text:
module SeaMass

import HDF5
import LibExpat

export MzmlbSpectrum, SmbSpectrum

type MzmlbSpectrum
  mzs
  intensities
end

function MzmlbSpectrum(filename::AbstractString, spectrumID::Number=1)
  HDF5.h5open(filename, "r") do file
    spectrumIndex = file["mzML_spectrumIndex"][spectrumID:spectrumID+1] + 1

    mzML = LibExpat.xp_parse(String(file["mzML"][spectrumIndex[1]:spectrumIndex[2]-1]))

    arrayLength = parse(Int, mzML["@defaultArrayLength"][1])
    if (arrayLength > 0)
      mzsDataset = mzML["binaryDataArrayList/binaryDataArray/cvParam[@accession='MS:1000514']/../binary"]["@externalDataset"][1]
      mzsOffset = parse(Int, mzML["binaryDataArrayList/binaryDataArray/cvParam[@accession='MS:1000514']/../binary"]["@offset"][1]) + 1
      intensitiesDataset = mzML["binaryDataArrayList/binaryDataArray/cvParam[@accession='MS:1000515']/../binary"]["@externalDataset"][1]
      intensitiesOffset = parse(Int, mzML["binaryDataArrayList/binaryDataArray/cvParam[@accession='MS:1000515']/../binary"]["@offset"][1]) + 1

      MzmlbSpectrum(
        file[mzsDataset][mzsOffset:mzsOffset+arrayLength-1],
        file[intensitiesDataset][intensitiesOffset:intensitiesOffset+arrayLength-1]
      )
    else
      MzmlbSpectrum(Float64[], Float32[])
    end
  end
end

type SmbSpectrum
  locations
  counts
  exposure
end

function SmbSpectrum(filename::AbstractString, spectrumID::Number=1)
  HDF5.h5open(filename, "r") do file
    if (HDF5.has(file,"countsIndex"))
      countsIndex = file["countsIndex"][spectrumID:spectrumID+1] + 1
    else
      countsIndex = [1 size(file["counts"])[1]+1]
    end

    if (countsIndex[2] - countsIndex[1] > 0)
      if (HDF5.has(file,"binLocations"))
        locations = file["binLocations"][countsIndex[1]+spectrumID-1:countsIndex[2]+spectrumID-1]
      else
        if(HDF5.has(file,"sampleLocations"))
          locations = file["sampleLocations"][countsIndex[1]:countsIndex[2]-1]
        else
          locations = file["centroidLocations"][countsIndex[1]:countsIndex[2]-1]
        end
      end

      if (HDF5.has(file,"exposures"))
        exposures = file["exposures"][spectrumID][1]
      else
        exposures = 1.0
      end

      SmbSpectrum(
        locations,
        file["counts"][countsIndex[1]:countsIndex[2]-1],
        exposures
      )
    else
       SmbSpectrum(Float64[], Float32[], Float32[])
    end
  end
end

end # module


text:
## Kernels
## =======

"Real+ -> [0, 1]"
kf1(x, β = 0.0001) = x / (x + β)
kf1β(β) = d -> kf1(d, β)

"(Log) Squared exponential kernel `α = 1/2l^2`, higher α is lower temperature  "
kse(d, α = 1000.0) = - α * d
kseα(α) = d -> kse(d, α)

"γ exponential"
γe(x, l = 1, γ = 1) = exp(-(x/l)^γ)

"Rational Quadratic"
rquad(x, l = 1, α = 1) = (1 + (x^2)/(2*α*l^2))^(-α)

"Power law kernel"
kpow(d, α = 1.0, k = 2) = -k * log(d) + log(α)

burr(x, c = 1, k = 40) =  log(c) + log(k) +  (c - 1) * log(x) - (k + 1)*log(1 + x^c)

# Temperature Modulation
const GLOBALKERNEL_ = Function[kse]

"Global Kernel"
function globalkernel!(k)
  global GLOBALKERNEL_
  GLOBALKERNEL_[1] = k
end

"Retrieve global kernel"
function globalkernel()
  global GLOBALKERNEL_
  GLOBALKERNEL_[1]
end

"Temporarily set global kernel"
function withkernel(thunk, k)
  globalkernel!(k)
  res = thunk()
  globalkernel!(kse)
  res
end

# (Cassette-based) Temperature Modulation
Cassette.@context AlphaCtx

"""
`f(args)` where temperature controlled with temperature `α`

```julia
x = normal(0.0, 1.0)
atα(10, rand, x ==ₛ 0.3)
```
"""
function atα(α, f, args...)
  ctx = AlphaCtx(metadata = α)
  Cassette.overdub(ctx, f, args...)
end

Cassette.overdub(ctx::AlphaCtx, ::typeof(kse), x, α) = kse(x, ctx.metadata)

"""

```julia
x = normal(0.0, 1.0)
@atα 100 rand(y ==ₛ 0.0)
```
"""
macro atα(a, fexpr)
  @pre fexpr.head == :call
  :(atα($(esc(a)), $(esc.(fexpr.args)...)))
end


text:
include("C:\\Users\\RR\\.julia\\dev\\Tokens\\test\\testsearch.jl")
using Profile
Profile.init(n=10^6,delay=0.0001)
@profile benchmark(128,300,10000000,"R:\\profile.csv")
open("r:\\prof.txt", "w") do s
    Profile.print(IOContext(s, :displaysize => (24, 500)))
end


text:
# (21 Jul 2017) Pccf Journal version
#  - 21 Jul 2017: keyword argumants in the cusum() function

function cusumStat(x, k)
    # cusum statistic
    n = length(x)
    stat = zeros(n)
    for i = 2:n
        stat[i] =  max(0, stat[i-1] + x[i] - k)
    end
    stat
end
# function cusum(x; update_width=10, expected_change=1.0, threshold=0.5)

function cusum(x; update_width=10, expected_change=1.0, threshold=0.5, tstoutput = false)
    n = length(x)
    stat_up = zeros(n)
    stat_down = zeros(n)
    wait_flag = 1
    chps = [1]
    update_points = []
    k_up, k_down = 0, 0
    mu0, mu1, mu2 = 0.0, 0.0, 0.0
    for i = 1 : n
        function update_params()
            train_dat = x[(i - update_width + 1) : i]
            mu0       = mean(train_dat)
            mu1       = mu0 + expected_change
            mu2       = mu0 - expected_change
            k_up      = (mu0 + mu1)/2
            k_down    = (mu0 + mu2)/2
        end
        if i - chps[end] == update_width && wait_flag == 1
            push!(update_points, i)
            update_params()
            wait_flag = 0
        end
        if wait_flag == 0
            stat_up[i]   = max(0, stat_up[i-1]   + x[i] - k_up)
            stat_down[i] = min(0, stat_down[i-1] + x[i] - k_down)
            # Change detection
            if (stat_up[i] >= threshold) || (abs(stat_down[i]) >= threshold)
                push!(chps, i)
                wait_flag = 1
            end
        end

    end
    if tstoutput
        return(Dict("chps" => chps, "update_points" => update_points, "stat_up" => stat_up, "stat_down" => stat_down))
    else
        return chps[2:end]
    end
end
#cusum(randn(100), 10, 1.0, 0.5)

function cusum_pccf(x, update_width, expected_change, threshold_min, threshold_max, mu, sigma, pccf_thr, tstoutput = false)
    n = length(x)
    probs = pccf_scale(Pccf.pccf(n, mu, sigma))
    hs = pccf_dynamic_h(probs, pccf_thr, threshold_min, threshold_max, false)
    chps = [1]

    stat_up = zeros(n)
    stat_down = zeros(n)
    wait_flag = 1
    update_points = []
    k_up, k_down = 0, 0
    mu0, mu1, mu2 = 0.0, 0.0, 0.0
    for i = 1 : n
        function update_params()
            train_dat = x[(i - update_width + 1) : i]
            mu0       = mean(train_dat)
            mu1       = mu0 + expected_change
            mu2       = mu0 - expected_change
            k_up      = (mu0 + mu1)/2
            k_down    = (mu0 + mu2)/2
        end
        if i - chps[end] == update_width && wait_flag == 1
            push!(update_points, i)
            update_params()
            wait_flag = 0
        end

        if wait_flag == 0
            stat_up[i]   = max(0, stat_up[i-1]   + x[i] - k_up)
            stat_down[i] = min(0, stat_down[i-1] + x[i] - k_down)
            # Change detection
            if (stat_up[i] > hs[i]) || (abs(stat_down[i]) > hs[i])
                push!(chps, i)
                wait_flag = 1
            end
        end

    end
    if tstoutput
        return(Dict("chps" => chps, "update_points" => update_points, "stat_up" => stat_up, "stat_down" => stat_down))
    else
        return chps[2:end]
    end
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
using Random, Test
using ModelParams

# Make test entryProb_jlcM with some zero entries
function make_test_entry_probs(rng, J, nc, nl)
    entryProb_jlcM = zeros(J, nl, nc);
    entryProb_jlM = rand(rng, J, nl);
    for j = 1 : J
        for l = 1 : nl
            prob_cV = rand(rng, nc);
            prob_cV[rand(rng, 1 : nc)] = 0.0;
            entryProb_jlcM[j, l, :] .= prob_cV .*
                (entryProb_jlM[j, l] / sum(prob_cV));
        end
    end
    return entryProb_jlcM
end

function test_value_cl(rng, nc, nl)
    return 1.0 .+ 2.0 .* rand(rng, nc, nl);
end

# Make type masses, setting some to 0
function type_masses_for_test(rng, J, nl)
    if nl == 1
        return 0.5 .+ rand(rng, J, nl);
    else
        typeMass_jlM = zeros(J, nl);
        for j = 1 : J
            typeMass_jlM[j, :] = 0.5 .+ rand(rng, nl);
            l = rand(rng, 1 : nl);
            typeMass_jlM[j, l] = 0.0;
        end
        return typeMass_jlM
    end
end

# College capacity matrix with given total capacity
# Some colleges do not occur in some places
function college_capacities_for_test(rng, nc, nl, totalCapacity)
    capacity_clM = rand(rng, nc, nl);
    if nl > 1
        for ic = 1 : nc
            l = rand(rng, 1 : nl);
            capacity_clM[ic, l] = 0.0;
        end
    end
    cSum  = sum(capacity_clM);
    capacity_clM .*= (totalCapacity / cSum);
    @assert isapprox(sum(capacity_clM), totalCapacity)
    return capacity_clM
end

function show_matrix(m :: Matrix{F1}; header = "Matrix:") where F1
    println(header);
    nr = size(m, 1);
    for ir = 1 : nr
        println(round.(m[ir,:], digits = 3));
    end
end

# Low capacity sequential ensures that there are full colleges
function test_entry_switches(J, nc; capacities = (:high, :low))
    if any(capacities .== :low)
        cap = 0.3 * J * nc;
    else
        cap = 100.0 * J * nc;
    end
    return [
        make_test_entry_sequ_multiloc(J, nc, nc + 1, cap),
        make_test_entry_sequ_multiloc(J, nc, nc + 1, cap; localOnlyV = [1]),
        make_entry_switches_oneloc(J, nc)
    ];
end

# Inputs are
# - `totalCapacity` of all colleges; in multiples of  total `typeMass`
function make_test_entry_sequ_multiloc(J, nc, nl, totalCapacity;
    localOnlyV = Vector{Int}())

    rng = MersenneTwister(45);
    typeMass_jlM = type_masses_for_test(rng, J, nl);
    totalMass = sum(typeMass_jlM);
    totalCapacity = 0.4 * totalMass;
    capacity_clM = college_capacities_for_test(rng, nc, nl, totalCapacity);
    if nl > 1
        valueLocal = 0.7;
        calValueLocal = true;
    else
        valueLocal = 0.0;
        calValueLocal = false;
    end

    switches = EntryDecisionSwitches(
        nTypes = J, nColleges = nc, nLocations = nl, valueLocal = valueLocal,
        typeMass_jlM = typeMass_jlM, capacity_clM = capacity_clM,
        entryPrefScale = 1.2,
        calValueLocal = calValueLocal);
    set_local_only_colleges!(switches, localOnlyV);
    @assert validate_es(switches)
    return switches
end

# ---------------------


text:
using TimeZones
using BenchmarkTools

suite = BenchmarkGroup()
suite["parse"] = BenchmarkGroup()

# Parsing a ZonedDateTime using the Z slot
function multiple(n)
    df = Dates.DateFormat("yyyymmddHH:MM:SS ZZZ")
    arr = Array{ZonedDateTime}(n)
    for (i, s) in enumerate(repeated("2016060701:02:03 America/Toronto", n))
        arr[i] = ZonedDateTime(s,df)
    end
    return arr
end
suite["parse"]["multiple"] = @benchmarkable multiple(1000)

suite


text:
using LinearAlgebra
A = [5 -2 -2; -2 4 -1; -2 -1 3]
L = cholesky(A)
##
println("A = ")
display(A)
println("L = ")
display(L.L)
##
@show norm(A - L.L*L.L')


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type DoubleValue
    value::Float64
    DoubleValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type DoubleValue
hash(v::DoubleValue) = ProtoBuf.protohash(v)
isequal(v1::DoubleValue, v2::DoubleValue) = ProtoBuf.protoisequal(v1, v2)
==(v1::DoubleValue, v2::DoubleValue) = ProtoBuf.protoeq(v1, v2)

type FloatValue
    value::Float32
    FloatValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type FloatValue
hash(v::FloatValue) = ProtoBuf.protohash(v)
isequal(v1::FloatValue, v2::FloatValue) = ProtoBuf.protoisequal(v1, v2)
==(v1::FloatValue, v2::FloatValue) = ProtoBuf.protoeq(v1, v2)

type Int64Value
    value::Int64
    Int64Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Int64Value
hash(v::Int64Value) = ProtoBuf.protohash(v)
isequal(v1::Int64Value, v2::Int64Value) = ProtoBuf.protoisequal(v1, v2)
==(v1::Int64Value, v2::Int64Value) = ProtoBuf.protoeq(v1, v2)

type UInt64Value
    value::UInt64
    UInt64Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type UInt64Value
hash(v::UInt64Value) = ProtoBuf.protohash(v)
isequal(v1::UInt64Value, v2::UInt64Value) = ProtoBuf.protoisequal(v1, v2)
==(v1::UInt64Value, v2::UInt64Value) = ProtoBuf.protoeq(v1, v2)

type Int32Value
    value::Int32
    Int32Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Int32Value
hash(v::Int32Value) = ProtoBuf.protohash(v)
isequal(v1::Int32Value, v2::Int32Value) = ProtoBuf.protoisequal(v1, v2)
==(v1::Int32Value, v2::Int32Value) = ProtoBuf.protoeq(v1, v2)

type UInt32Value
    value::UInt32
    UInt32Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type UInt32Value
hash(v::UInt32Value) = ProtoBuf.protohash(v)
isequal(v1::UInt32Value, v2::UInt32Value) = ProtoBuf.protoisequal(v1, v2)
==(v1::UInt32Value, v2::UInt32Value) = ProtoBuf.protoeq(v1, v2)

type BoolValue
    value::Bool
    BoolValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type BoolValue
hash(v::BoolValue) = ProtoBuf.protohash(v)
isequal(v1::BoolValue, v2::BoolValue) = ProtoBuf.protoisequal(v1, v2)
==(v1::BoolValue, v2::BoolValue) = ProtoBuf.protoeq(v1, v2)

type StringValue
    value::AbstractString
    StringValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type StringValue
hash(v::StringValue) = ProtoBuf.protohash(v)
isequal(v1::StringValue, v2::StringValue) = ProtoBuf.protoisequal(v1, v2)
==(v1::StringValue, v2::StringValue) = ProtoBuf.protoeq(v1, v2)

type BytesValue
    value::Array{UInt8,1}
    BytesValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type BytesValue
hash(v::BytesValue) = ProtoBuf.protohash(v)
isequal(v1::BytesValue, v2::BytesValue) = ProtoBuf.protoisequal(v1, v2)
==(v1::BytesValue, v2::BytesValue) = ProtoBuf.protoeq(v1, v2)

export DoubleValue, FloatValue, Int64Value, UInt64Value, Int32Value, UInt32Value, BoolValue, StringValue, BytesValue


text:
############################################################
# Affinity Propagation
############################################################

type AffinityPropagationOpts
    max_iter     ::Int
    n_stop_check ::Int
    damp         ::Float64
    display      ::Symbol
end

function AffinityPropagationOpts(;
    max_iter::Integer = 500,    # max number of iterations
    n_stop_check::Integer = 10, # stop if exemplars not changed for this number of iterations
    damp::FloatingPoint = 0.5,  # damping factor for message updating, 0 means no damping
    display::Symbol = :iter     # verboseness: :iter, :final, :none
    )

    AffinityPropagationOpts(int(max_iter), int(n_stop_check), float(damp), display)
end

type AffinityPropagationResult
    exemplar_index ::Vector{Int} # index for exemplars (centers)
    assignments    ::Vector{Int} # assignments for each point
    iterations     ::Int         # number of iterations executed
    converged      ::Bool        # converged or not
end

function affinity_propagation{T<:FloatingPoint}(S::Matrix{T}, opts::AffinityPropagationOpts)
    n = size(S, 1)
    converged = false

    exemplars = falses(n, opts.n_stop_check)
    psi = zeros(eltype(S), n, n)
    phi = zeros(eltype(S), n, n)

    n_iter = 0
    IC = nothing

    for n_iter = 1:opts.max_iter
        affprop_update_message!(psi, phi, S, n, opts)

        IC = diag(psi) + diag(phi) .> 0
        exemplars[:, mod(n_iter-1, opts.n_stop_check)+1] = IC

        if opts.display == :iter
            @printf("%7d: %3d exemplars identified\n", n_iter, nnz(IC))
        end

        if n_iter > opts.n_stop_check
            converged = true
            for i = 1:opts.n_stop_check
                if IC != exemplars[:, i]
                    converged = false
                    break
                end
            end
            if converged
                break
            end
        end
    end

    assignment, exemplar_index = affprop_decide_assignment(phi, S, n, IC)
    if opts.display == :iter || opts.display == :final
        if converged
            println("affinity propagation converged with $n_iter iterations")
        else
            println("affinity propagation terminated without convergence after $n_iter iterations")
        end
    end

    return AffinityPropagationResult(exemplar_index, assignment, n_iter, converged)
end

affinity_propagation(x::Matrix; opts...) = affinity_propagation(x, AffinityPropagationOpts(;opts...))

############################################################
# Message Updating
#
# Loopy message passing sometimes can be tricky. Parallel
# updating vs. in-place updating; damping vs. no-damping;
# as well as the order of message updating can affect the
# convergence behavior drastically.
############################################################
function affprop_update_message!{T<:FloatingPoint}(psi::Matrix{T}, phi::Matrix{T}, S::Matrix{T}, n, opts::AffinityPropagationOpts)
    # update phi
    for i = 1:n
        SM = S[i,:] + psi[i,:]
        I = indmax(SM)
        Y = SM[I]
        SM[I] = -Inf
        Y2 = max(SM)
        val = S[i,:] - Y;
        val[I] = S[i,I] - Y2;
        phi[i,:] = opts.damp*phi[i,:] + (1-opts.damp)*val
    end

    # update psi
    for j = 1:n
        RP = phi[:,j]
        idx = 1:n .!= j
        RP[idx] = max(RP[idx], 0)
        val = sum(RP) - RP;
        val[idx] = min(val[idx], 0)
        psi[:,j] = opts.damp*psi[:,j] + (1 - opts.damp)*val
    end
end

function affprop_decide_assignment{T<:FloatingPoint}(phi::Matrix{T}, S::Matrix{T}, n::Integer, exemplar::BitArray{1})
    assignment = zeros(Int, n)
    Iexp = find(exemplar)

    for i = 1:n
        I = indmax(phi[i,Iexp] + S[i,Iexp])
        assignment[i] = Iexp[I]
    end

    # make sure exemplars are assigned to themselves
    assignment[Iexp] = Iexp

    return assignment, Iexp
end


text:
using TimeZones
function read_and_capitalize(f::IOStream)
    ss = uppercase(read(f, String))
    println(ss)
    return ss
end
f = open("hello.txt","w")
write(f,"Hi, MarsBighead")
close(f)
f = open("hello.txt")
read_and_capitalize(f)
close(f)
d = DateTime(2018)
println(d)

max_stars_repo_name: MarsBighead/mustang

max_issues_repo_name: MarsBighead/mustang

max_forks_repo_name: MarsBighead/mustang


text:
############################################################################################
"""
$(TYPEDEF)
Metropolis algorithm container.

# Fields
$(TYPEDFIELDS)
"""
mutable struct Metropolis{R<:ℓDensityResult} <: MCMCKernel
    "Cached Result of last propagation step."
    result::R
    function Metropolis(result::R) where {R<:ℓObjectiveResult}
        return new{R}(result)
    end
end
function Metropolis(sym::S; kwargs...) where {S<:Union{Symbol,NTuple{k,Symbol} where k}}
    return MCMC(Metropolis, sym; kwargs...)
end

function update!(kernel::Metropolis, objective::Objective, up::BaytesCore.UpdateTrue)
    ## Update log-target result with current (latent) data
    kernel.result = ModelWrappers.ℓDensityResult(objective)
    return nothing
end
function update!(kernel::Metropolis, objective::Objective, up::BaytesCore.UpdateFalse)
    return nothing
end
############################################################################################
"""
$(TYPEDEF)
Representation of a trajectory.

# Fields
$(TYPEDFIELDS)
"""
struct TrajectoryMetropolis{
    S<:AbstractMatrix,F<:ModelWrappers.ℓObjectiveResult,T<:AbstractFloat
}
    "Proposal Covariance"
    Σ::S
    "Log Target result at initial point, see ModelWrappersesTools."
    result₀::F
    "Discretization size"
    ϵ::T
    function TrajectoryMetropolis(
        Σ::S, result₀::F, ϵ::T
    ) where {S<:AbstractMatrix,F<:ModelWrappers.ℓObjectiveResult,T<:AbstractFloat}
        return new{S,F,T}(Σ, result₀, ϵ)
    end
end

function move(_rng::Random.AbstractRNG, trajectory::T) where {T<:TrajectoryMetropolis}
    @unpack result₀, Σ, ϵ = trajectory
    return rand(_rng, MvNormal(result₀.θᵤ, ϵ .* Σ))
end

function checkfinite(
    trajectory::TrajectoryMetropolis, result::ModelWrappers.ℓObjectiveResult
)
    return checkfinite(trajectory.result₀.ℓθᵤ, result.ℓθᵤ, result)
end

############################################################################################
"propagate Metropolis sampler forward"
function propagate(
    _rng::Random.AbstractRNG, kernel::Metropolis, tune::MCMCTune, objective::Objective
)
    @unpack result = kernel
    @unpack Σ = tune.proposal
    @unpack ϵ = tune.stepsize
    ## Create new trajectory
    trajectory = TrajectoryMetropolis(Σ, kernel.result, ϵ)
    ## Make Proposal step
    resultᵖ = ModelWrappers.log_density(objective, move(_rng, trajectory))
    ## Check if proposal diverges
    divergent = !checkfinite(result, resultᵖ)
    if divergent
        return resultᵖ,
        divergent, BaytesCore.AcceptStatistic(zero(typeof(ϵ)), false),
        DiagnosticsMetropolis(ϵ)
    end
    ## Calculate Proposal density and acceptance rate ~ not needed for symmetric proposals
    #ℓqᵤ = logpdf(MvNormal(resultᵖ.θᵤ, ϵ * Σ), result.θᵤ)
    #ℓqᵤᵖ = logpdf(MvNormal(result.θᵤ, ϵ * Σ), resultᵖ.θᵤ)
    ## Pack and return output
    accept_statistic = BaytesCore.AcceptStatistic(_rng, (resultᵖ.ℓθᵤ - result.ℓθᵤ)) #+ (ℓqᵤ - ℓqᵤᵖ))
    return resultᵖ, divergent, accept_statistic, DiagnosticsMetropolis(ϵ)
end

function get_acceptrate(
    _rng::Random.AbstractRNG, kernel::Metropolis, objective::Objective, Σ::M
) where {M<:AbstractMatrix}
    ## Current logposterior and proposal density
    @unpack result = kernel
    ## Function of stepsize
    return function acceptrate(ϵ::T) where {T<:Real}
        trajectory = TrajectoryMetropolis(Σ, result, ϵ)
        ## Make Proposal step
        resultᵖ = ModelWrappers.log_density(objective, move(_rng, trajectory))
        #ℓqᵤ = logpdf(MvNormal(resultᵖ.θᵤ, ϵ * Σ), result.θᵤ)
        #ℓqᵤᵖ = logpdf(MvNormal(result.θᵤ, ϵ * Σ), resultᵖ.θᵤ)
        ## Return acceptance rate
        return exp((resultᵖ.ℓθᵤ - result.ℓθᵤ)) #+ (ℓqᵤ - ℓqᵤᵖ)) #Unbounded accpetance rate
    end
end

############################################################################################
# Export
export Metropolis


text:
module BenchmarkingBrentFmin

using Printf
using BenchmarkTools

using OptimPackNextGen.Brent
using Brent: fmin_atol, fmin_rtol

include("fmin-funcs.jl")
using .FminTestFunctions:
    TestFunc,
    incr_nevals,
    reset_nevals,
    get_nevals,
    func,
    ident,
    local_lower,
    local_upper,
    global_lower,
    global_upper,
    xsol,
    fsol,
    ampgo_2,
    ampgo_3,
    ampgo_4,
    ampgo_5,
    ampgo_6,
    ampgo_7,
    ampgo_8,
    ampgo_9,
    ampgo_10,
    ampgo_11,
    ampgo_12,
    ampgo_13,
    ampgo_14,
    ampgo_15,
    ampgo_18,
    ampgo_20,
    ampgo_21,
    ampgo_22,
    gsl_fmin_1,
    gsl_fmin_2,
    gsl_fmin_3,
    gsl_fmin_4

function run(;
             T::Type{<:AbstractFloat}=Float64,
             prec::Bool=false, # tests with high precision
             )
    if prec
        atol = 1e-100
        rtol = 1e-15
    else
        atol = fmin_atol(T)
        rtol = fmin_rtol(T)
    end
    for _f in (ampgo_2, ampgo_3, ampgo_4, ampgo_5, ampgo_6, ampgo_7,
               ampgo_8, ampgo_9, ampgo_10, ampgo_11, ampgo_12, ampgo_13,
               ampgo_14, ampgo_15, ampgo_18, ampgo_20, ampgo_21, ampgo_22,
               gsl_fmin_1, gsl_fmin_2, gsl_fmin_3, gsl_fmin_4)
        f = TestFunc{T}(_f)
        a = local_lower(f)
        b = local_upper(f)
        x0 = xsol(f)
        f0 = fsol(f)
        reset_nevals()
        (xm, fm, lo, hi) = fmin(f, a, b; rtol=rtol, atol=atol)
        nevals = get_nevals()
        println("Problem $(ident(f)) (T=$T):")
        @printf(" ├─ nevals = %d\n", nevals)
        if prec
            @printf(" ├─ x0 = %.13f, f0 = %.13f\n", xm, fm)
            @printf(" └─ prec = %.3e\n", (hi - lo)/2)
        else
            @printf(" ├─")
            @btime fmin($f, $a, $b);
            @printf(" ├─ xm     = %.15f ± %.3e\n", xm, (hi - lo)/2)
            @printf(" ├─ x0     = %.15f\n", x0)
            @printf(" ├─ x_err  = %.3e\n", abs(xm - x0))
            @printf(" ├─ f(xm)  = %.15f\n", fm)
            @printf(" ├─ f(x0)  = %.15f\n", f0)
            @printf(" └─ f_err  = %.3e\n", abs(fm - f0))
        end
    end
end

end # module


text:
const TTMCS = TMazeCumulantSchedules

mutable struct Cumulants <: LoggerKeyData
    cs::Array{Float64,2}
    agent_pos::Array{Float64,2}
    function Cumulants(logger_init_info)
        num_logged_steps = logger_init_info[LoggerInitKey.TOTAL_STEPS]
        new(zeros(4,num_logged_steps),zeros(2,num_logged_steps))
    end
end

function lg_step!(self::Cumulants, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    true_values = TTMCS.get_cumulant_eval_values(env.cumulant_schedule)
    if s_next[3] != 0.0
        true_values[1] = s_next[3]
    end
    true_values[3] = s_next[5]
    self.cs[:,cur_step_total] = true_values
    self.agent_pos[:,cur_step_total] = s_next[1:2]
end

function lg_episode_end!(self::Cumulants, cur_step_in_episode, cur_step_total)
end

function save_log(self::Cumulants, save_dict::Dict)
    save_dict[:cumulants] = self.cs
    save_dict[:agent_pos] = self.agent_pos
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using Test
include("../src/SymbolicMajoranas.jl")

@testset "BareCoefficient" begin

    # for convenience
    BC = BareCoefficient

    # check order defined on bare coefficients
    @test BC(-3, "a") < BC(2, "d")
    @test BC(2, "b") < BC(2, "d")
    @test BC(2, "b") == BC(2, "b")

end

@testset "Coefficient" begin

    # for convenience
    BC = BareCoefficient
    C = Coefficient

    # check in place multiplication by a numerical factor
    @test begin

        c = C(-2//5, [BC(1, "b"), BC(2, "c")], [BC(0, "a")])
        times!(c, -1//2)
        c.num == 1//5

    end

    # check zeroing a coefficent works as intended
    @test begin

        c = C(-2//5, [BC(1, "b"), BC(2, "c")], [BC(0, "a")])
        zero!(c)
        c.num == 0 && isempty(c.top) && isempty(c.bot)

    end

    # check canonicalization and addition of coefficients
    @test begin

        c1 = C(1, [BC(1, "b"), BC(2, "c"), BC(4, "e")], [BC(4, "e"), BC(0, "a")])
        c2 = C(2, [BC(2, "c"), BC(1, "b")], [BC(0, "a")])
        canonicalize!(c1)
        canonicalize!(c2)
        addable(c1, c2)
        c3 = c1 + c2
        c3.num == 3

    end

end

@testset "CoefficientSum" begin

    # for convenience
    BC = BareCoefficient
    C = Coefficient
    CS = CoefficientSum

    # check the simplify! method
    @test begin

        cs = CS([
            C(1, [BC(1, "b"), BC(2, "c")], [BC(0, "a")]),
            C(2, [BC(1, "b"), BC(2, "c")], [BC(0, "a")]),
            C(3, [BC(3, "d")], []),
            C(4, [BC(3, "d"), BC(0, "z")], [BC(0, "z")]),
            zero(C),
        ])

        simplify!(cs)
        size(cs) == (2,) && cs[1].num == 3 && cs[2].num == 7

    end

    # check in place multiplication by a scalar
    @test begin

        cs = CS([
            C(1, [BC(1, "b"), BC(2, "c")], [BC(0, "a")]),
            C(4, [BC(3, "d"), BC(0, "z")], [BC(0, "z")]),
        ])

        times!.(cs, -1im)
        cs[1].num == -1im && cs[2].num == -4im

    end

end

@testset "MajoranaProduct" begin

    # for convenience
    MP = MajoranaProduct

    # test canonicalization
    @test begin

        mp = MP([3, 1, 2, 1])
        ic1 = iscanonical(mp)
        n = canonicalize!(mp)
        ic2 = iscanonical(mp)

        !ic1 && n == 1 && ic2 && length(mp) == 2

    end

    # test * and ==
    @test begin

        mp1 = MP([1, 2, 3, 4])
        mp2 = MP([3, 1])
        mp3, n = mp1 * mp2
        mp3 == MP([2, 4]) && n == -1

    end

    # test commute
    @test commute(MP([1, 2]), MP([2, 4])) == false
    @test commute(MP([1, 2]), MP([3, 4])) == true

end

@testset "Term" begin

    # for convenience
    BC = BareCoefficient
    C = Coefficient
    CS = CoefficientSum
    MP = MajoranaProduct
    T = Term

    # TODO add tests here

end

@testset "TermSum" begin

    # for convenience
    BC = BareCoefficient
    C = Coefficient
    CS = CoefficientSum
    MP = MajoranaProduct
    T = Term
    TS = TermSum

    # TODO add tests here

end

max_stars_repo_name: aormorningstar/SymbolicMajoranas

max_issues_repo_name: aormorningstar/SymbolicMajoranas

max_forks_repo_name: aormorningstar/SymbolicMajoranas


text:
export UCBExplorer

using Random
using Flux

Base.@kwdef mutable struct UCBExplorer{R<:AbstractRNG} <: AbstractExplorer
    c::Float64
    actioncounts::Vector{Float64}
    step::Int
    rng::R
    is_training::Bool = true
end

"""
    UCBExplorer(na; c=2.0, ϵ=1e-10, step=1, seed=nothing)

# Arguments
- `na` is the number of actions used to create a internal counter.
- `t` is used to store current time step.
- `c` is used to control the degree of exploration.
- `seed`, set the seed of inner RNG.
- `is_training=true`, in training mode, time step and counter will not be updated.
"""
UCBExplorer(na; c = 2.0, ϵ = 1e-10, step = 1, rng = Random.GLOBAL_RNG, is_training = true) =
    UCBExplorer(c, fill(ϵ, na), 1, rng, is_training)

function (p::UCBExplorer)(values::AbstractArray)
    v, inds = find_all_max(@. values + p.c * sqrt(log(p.step + 1) / p.actioncounts))
    action = sample(p.rng, inds)
    if p.is_training
        p.actioncounts[action] += 1
        p.step += 1
    end
    action
end


text:
module PandaModels

import JuMP

import InfrastructureModels; const _IM = InfrastructureModels

import PowerModels; const _PM = PowerModels
import PowerModels: ids, ref, var, con, sol, nw_id_default

import Memento
const LOGGER = Memento.getlogger(PowerModels)

import JSON

import Cbc
import Ipopt
import Juniper

export run_powermodels_opf,
    run_powermodels_ots,
    run_powermodels_tnep,
    run_powermodels_powerflow,
    run_powermodels_mn_storage,
    run_powermodels_custom

include("input/pp_to_pm.jl")
# include("input/time_series.jl")
# include("input/tools.jl")
include("models/call_powermodels.jl")

end


text:
function _write_coverage(filename::Symbol;
                         dump_coverage::Bool = false,
                         dump_coverage_io::IO = stdout,
                         min_padding::Integer,
                         max_padding::Integer)::Nothing
    coverage_filename = string(filename, ".microcov")
    coverage_lines = _compute_coverage_file_contents(filename;
                                                     min_padding = min_padding,
                                                     max_padding = max_padding)
    @debug("Writing coverage to file",
           source_file = filename,
           coverage_file = coverage_filename)
    rm(coverage_filename; force = true, recursive = true)
    open(coverage_filename, "w") do io
        if dump_coverage
            # println(dump_coverage_io, "\n")
            println(dump_coverage_io, "# $(coverage_filename):")
        end
        for x in coverage_lines
            if dump_coverage
                println(dump_coverage_io, x)
            end
            println(io, x)
        end
        if dump_coverage
            # println(dump_coverage_io, "\n")
        end
    end
    return nothing
end


text:
include("update_p.jl")
include("update_beta.jl")
include("update_gamma.jl")
include("update_eta.jl")
include("update_lambda.jl")
include("update_mu.jl")
include("update_nu.jl")
include("update_omega.jl")
include("update_tau.jl")
include("update_psi.jl")
include("update_v.jl")
include("update_zeta.jl")
include("update_sigma.jl")
include("update_phi.jl")

function update_state!(state::State, data::Data, prior::Prior, tuners::Tuners;
                       fix::Vector{Symbol}=Symbol[],
                       flags::Vector{Flag}=Flag[])
  # Returns true if a parameter is specified as fixed (to not update).
  isfixed(sym::Symbol) = sym in fix

  # NOTE: Update parameters if they are not held fixed, as specified in `fix`.
  isfixed(:p) || update_p!(state, data, prior)
  isfixed(:beta) || update_beta!(state, data, prior, flags)

  # NOTE: `lambda` must be updated immediately after updating beta since lambda
  # is marginalized over in the update of beta. Marginalizing over lambda makes
  # updates for beta slower, but the mixing will be faster.
  isfixed(:lambdaC) || isfixed(:lambda) || update_lambdaC!(state, data, prior, flags)
  isfixed(:lambdaT) || isfixed(:lambda) || update_lambdaT!(state, data, prior, flags)

  # NOTE: If v and zeta are marginalized over in the updates of lambda,
  # then they must be updated immediately after updating lambda.
  isfixed(:vC) || isfixed(:v) || update_v!('C', state, data, prior)
  isfixed(:vT) || isfixed(:v) || update_v!('T', state, data, prior)
  isfixed(:zetaC) || isfixed(:zeta) || update_zeta!('C', state, data, prior)
  isfixed(:zetaT) || isfixed(:zeta) || update_zeta!('T', state, data, prior)

  isfixed(:gammaC) || isfixed(:gamma) || update_gammaC!(state, data, prior)
  isfixed(:gammaT) || isfixed(:gamma) || update_gammaT!(state, data, prior)
  isfixed(:etaC) || isfixed(:eta) ||  update_etaC!(state, data, prior)
  isfixed(:etaT) || isfixed(:eta) || update_etaT!(state, data, prior)

  isfixed(:mu) || update_mu!(state, data, prior)
  isfixed(:tau) || update_tau!(state, data, prior)
  isfixed(:omega) || update_omega!(state, data, prior)
  isfixed(:psi) || update_psi!(state, data, prior)
  isfixed(:nu) || update_nu!(state, data, prior, tuners)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
includet("aidpolicy.jl")
using Plots
#6 (or maybe 7) is the limit on my laptop
N=3
G=6
pi_t = 10
pi_b = 1
nsims = pi_t - pi_b+1
obj_opt = zeros(nsims)
obj_tic = zeros(nsims)

for (pi, i) in zip(pi_b:pi_t, 1:nsims)
    tree = initializeGame(N, G)
    calculate_TIC(tree, pi)
    obj_tic[i] = get_objective(tree)
    obj_opt[i] = optimalRule(tree, pi)
end

obj_unif = zeros(nsims) .+ (1- 1/(N-1))^2 .+ (N-2)*(0-1/(N-1))^2
obj_ante = zeros(nsims) .+ (1- 1/(N))^2 .+ (N-1)*(0-1/(N))^2
plot(pi_b:pi_t,[ obj_ante obj_unif obj_opt obj_tic], ylims=[0,1.1], xlabel="B", ylabel="Objective",label=["Ex-Ante Uniform" "Ex-Post Uniform" "Optimal" "T-IC"])

#weights plot
function weights_plot()
    N=3
    G=6
    pi=10
    tic_weights = zeros(G)
    tree = initializeGame(N, G)
    calculate_TIC(tree, pi)
    obj, opt_weights = optimalRule(tree, pi, true, true)
    obj, no_inc_weights = optimalRule(tree, pi, false, true)
    getWeights!(tree, tic_weights, 0)
    plot([ no_inc_weights opt_weights tic_weights], ylim = [0.08, 0.15], xlabel="Period", ylabel="Weight",labels=["Optimal, No IC Constraints" "Optimal Constrained" "T-IC "])
    savefig("figures/weights10.pdf")
end

weights_plot()

max_stars_repo_name: evanmunro/draft-policy

max_issues_repo_name: evanmunro/draft-policy

max_forks_repo_name: evanmunro/draft-policy


text:
using Test
using Logging
using Distributions: Normal, Poisson, Gamma, LogNormal, Bernoulli, Geometric, MvNormal, truncated, Dirichlet, Categorical, logpdf, MvLogNormal, DiscreteUniform, fit, Beta, Binomial, DiscreteNonParametric, rand
using StatsBase: mean, std
using PrettyPrint: pprintln
using Random
using PDMats
using JSON

using CrimsonSkyline

const Maybe{T} = Union{T,Nothing}

Random.seed!(2021)

include("basic.jl")
include("condition_and_observe.jl")
include("distributions.jl")
include("field.jl")
include("io.jl")
include("forward.jl")
include("importance_sampling.jl")
include("metropolis.jl")
include("nested.jl")
include("typed_trace.jl")
include("statistics.jl")
include("db.jl")
include("kernel.jl")
include("dsl.jl")


text:
#
# returns the indexes (i,j) of given the index of the pilled-up columns
#
function setindex2D(ipair,npairs)
  j = ipair%npairs
  if j == 0
    j = 1
  end
  i = round(Int64,(npairs-j)/npairs)+1
  return i, j
end

function setipair2D(n,i,j)
  if j >= i
    ipair = (i-1)*n + j
  else
    ipair = (j-1)*n + i
  end
  ipair = ipair - round(Int64,(i^2+i)/2)
  return ipair
end