text:
# This script can be used to simulate the onboard-Kalman filter using real recorded sensor data.
# This can be useful to find the optimal Kalman parameters (Q,R) or to just do experiments and find an optimal sensor fusion process.

using PyPlot
using JLD

Q = diagm([0.1,0.1,0.1,0.1,0.1])
R = diagm([0.1,0.1,0.1,10.0])

function main(code::AbstractString)
    global Q, R
    log_path_record = "$(homedir())/open_loop/output-record-$(code).jld"
    d_rec = load(log_path_record)

    imu_meas    = d_rec["imu_meas"]
    gps_meas    = d_rec["gps_meas"]
    cmd_log     = d_rec["cmd_log"]
    cmd_pwm_log = d_rec["cmd_pwm_log"]
    vel_est     = d_rec["vel_est"]
    pos_info    = d_rec["pos_info"]

    t0      = max(imu_meas.t[1],vel_est.t[1],gps_meas.t[1])
    t_end   = min(imu_meas.t[end],vel_est.t[end],gps_meas.t[end],cmd_pwm_log.t[end])

    l_A = 0.125
    l_B = 0.125

    dt = 1/25
    t = t0:dt:t_end
    sz = length(t)
    y = zeros(sz,4)
    u = zeros(sz,2)

    P = zeros(5,5)
    x_est = zeros(length(t),5)

    yaw0 = imu_meas.z[1,6]
    gps_dist = zeros(length(t))

    yaw_prev = yaw0

    for i=2:length(t)
        # Collect measurements and inputs for this iteration
        y_gps = gps_meas.z[gps_meas.t.>t[i],:][1,:]
        y_yaw = imu_meas.z[imu_meas.t.>t[i],6][1]-yaw0
        #a_x   = imu_meas.z[imu_meas.t.>t[i],7][1]
        y_vel_est = vel_est.z[vel_est.t.>t[i]][1]
        y[i,:] = [y_gps y_yaw y_vel_est]
        y[:,3] = unwrap!(y[:,3])
        u[i,:] = cmd_log.z[cmd_log.t.>t[i],:][1,:]
        #u[i,1] = (cmd_pwm_log.z[cmd_pwm_log.t.>t[i],1][1]-94.14)/2.7678
        #u[i,2] = (cmd_pwm_log.z[cmd_pwm_log.t.>t[i],2][1]-91.365)/105.6
        # if cmd_pwm_log.z[cmd_pwm_log.t.>t[i],1][1] == 90
        #     u[i,1] = 0
        # end
        # if cmd_pwm_log.z[cmd_pwm_log.t.>t[i],2][1] == 90
        #     u[i,2] = 0
        # end

        # Adapt R-value of GPS according to distance to last point:
        gps_dist[i] = norm(y[i,1:2]-x_est[i-1,1:2])
        if gps_dist[i] > 0.3
            R[1,1] = 1+100*gps_dist[i]^2
            R[2,2] = 1+100*gps_dist[i]^2
        else
            R[1,1] = 1
            R[2,2] = 1
        end

        args = (u[i,:],dt,l_A,l_B)

        # Calculate new estimate
        (x_est[i,:], P) = ekf(simModel,x_est[i-1,:]',P,h,y[i,:]',Q,R,args)
    end
    # ax1=subplot(211)
    # plot(t,y,"-x",t,x_est,"-*")
    # grid("on")
    # legend(["x","y","psi","v"])
    # subplot(212,sharex=ax1)
    # plot(t,u)
    # grid("on")
    # legend(["a_x","d_f"])
    figure(1)
    ax=subplot(5,1,1)
    for i=1:4
        subplot(5,1,i,sharex=ax)
        plot(t,y[:,i],t,x_est[:,i],"-*")
        grid("on")
    end
    subplot(5,1,5,sharex=ax)
    plot(cmd_pwm_log.t,cmd_pwm_log.z)
    grid("on")

    figure(2)
    subplot(2,1,1)
    title("Comparison raw GPS data and estimate")
    plot(t-t0,y[:,1],t-t0,x_est[:,1],"-*")
    xlabel("t [s]")
    ylabel("Position [m]")
    legend(["x_m","x_est"])
    grid("on")
    subplot(2,1,2)
    plot(t-t0,y[:,2],t-t0,x_est[:,2],"-*")
    xlabel("t [s]")
    ylabel("Position [m]")
    legend(["y_m","y_est"])
    grid("on")

    figure(3)
    plot(t,gps_dist)
    title("GPS distances")
    grid("on")
    # figure()
    # plot(gps_meas.z[:,1],gps_meas.z[:,2],"-x",x_est[:,1],x_est[:,2],"-*")
    # grid("on")
    # title("x-y-view")
    figure(4)
    plot(t-t0,x_est,"-x",pos_info.t-t0,pos_info.z[:,[6,7,4,10,16]],"-*")
    title("Comparison simulation (-x) onboard-estimate (-*)")
    grid("on")
    legend(["x","y","psi","v","psi_drift"])

    unwrapped_yaw_meas = unwrap!(imu_meas.z[:,6])
    figure(5)
    plot(t-t0,x_est[:,3],"-*",imu_meas.t-t0,unwrapped_yaw_meas-unwrapped_yaw_meas[1])
    title("Comparison yaw")
    grid("on")

    figure(6)
    plot(t-t0,x_est[:,4],"-*",vel_est.t-t0,vel_est.z)
    grid("on")
    title("Comparison of velocity measurement and estimate")
    legend(["estimate","measurement"])

    #figure(7)
    #plt[:hist](gps_dist,300)
    #grid("on")
    nothing
end

function h(x,args)
    C = [eye(4) zeros(4,1)]
    C[4,4] = 0
    #C[3,3] = 0
    C[3,5] = 1
    return C*x
end
function ekf(f, mx_k, P_k, h, y_kp1, Q, R, args)
    xDim    = size(mx_k,1)
    mx_kp1  = f(mx_k,args)
    A       = numerical_jac(f,mx_k,args)
    P_kp1   = A*P_k*A' + Q
    my_kp1  = h(mx_kp1,args)
    H       = numerical_jac(h,mx_kp1,args)
    P12     = P_kp1*H'
    K       = P12*inv(H*P12+R)
    mx_kp1  = mx_kp1 + K*(y_kp1-my_kp1)
    P_kp1   = (K*R*K' + (eye(xDim)-K*H)*P_kp1)*(eye(xDim)-K*H)'
    return mx_kp1, P_kp1
end

function simModel(z,args)
   # kinematic bicycle model
   # u[1] = acceleration
   # u[2] = steering angle
    (u,dt,l_A,l_B) = args
    bta = atan(l_A/(l_A+l_B)*tan(u[2]))

    zNext = copy(z)
    zNext[1] = z[1] + dt*(z[4]*cos(z[3] + bta))                     # x
    zNext[2] = z[2] + dt*(z[4]*sin(z[3] + bta))                     # y
    zNext[3] = z[3] + dt*(z[4]/l_B*sin(bta))                        # psi
    zNext[4] = z[4] + dt*(u[1] - 0.63 * z[4]^2 * sign(z[4]))       # v
    #zNext[4] = z[4] + dt*(z[6])                                     # v
    #zNext[4] = z[4] + dt*(sqrt(z[6]^2+z[7]^2))                     # v
    zNext[5] = z[5]                                                 # psi_drift
    return zNext
end

function numerical_jac(f,x,args)
    y = f(x,args)
    jac = zeros(size(y,1),size(x,1))
    eps = 1e-5
    xp = copy(x)
    for i = 1:size(x,1)
        xp[i] = x[i]+eps/2.0
        yhi = f(xp,args)
        xp[i] = x[i]-eps/2.0
        ylo = f(xp,args)
        xp[i] = x[i]
        jac[:,i] = (yhi-ylo)/eps
    end
    return jac
end

function initPlot()             # Initialize Plots for export
    linewidth = 0.4
    rc("axes", linewidth=linewidth)
    rc("lines", linewidth=linewidth, markersize=2)
    #rc("font", family="")
    rc("axes", titlesize="small", labelsize="small")        # can be set in Latex
    rc("xtick", labelsize="x-small")
    rc("xtick.major", width=linewidth/2)
    rc("ytick", labelsize="x-small")
    rc("ytick.major", width=linewidth/2)
    rc("legend", fontsize="small")
    rc("font",family="serif")
    rc("font",size=8)
    rc("figure",figsize=[4.5,3])
    #rc("pgf", texsystem="pdflatex",preamble=L"""\usepackage[utf8x]{inputenc}\usepackage[T1]{fontenc}\usepackage{lmodern}""")
end

function unwrap!(p)
    length(p) < 2 && return p
    for i = 2:length(p)
        d = p[i] - p[i-1]
        if abs(d) > pi
            p[i] -= floor((d+pi) / (2*pi)) * 2pi
        end
    end
    return p
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
# Minimal experiment intended to test the experiment runner
# Should plot x^2

using SubmodularMaximization
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics
using Printf
using Random

close("all")

experiment_name = "minimal_experiment"
data_folder = "./data"
reprocess = true

scales = 1:100
repetitions = 10

trial_fun(x) = mean(y->rand(x, x), 1:repetitions)
print_summary(x) = println("Square: ", x)

results = run_experiments(scales,
                          trial_fun=trial_fun,
                          print_summary=print_summary,
                          experiment_name=experiment_name,
                          data_folder=data_folder,
                          reprocess=reprocess
                         )

y = [sum(results[x]) for x in scales]

plot(scales, y)

ext: jl

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
"""
printhist(tablein,headers;TeX=false,figures=5)

You are not supposed to look at this code. It is mostly repulsive
bookkeeping.

Inputs:
tablein: columns of equal length with the data
         Pad the short ones with NaNs. My codes which all this
         do that. The formatting turns the NaNs into spaces of
         the correct length.
headers: the titles for the columns
         example: headers=["foo", "bar"]
         Do not add the header for the iteration counter. I do that.
TeX:     Format the tabel in LaTeX or not?
figures: The printf format will be (figures+7).figures e
         So if figures = 5 the data will be formatted in %12e5.
         This makes everything line up.

Compare iteration histories from SIAMFANLEquations family of solvers.
This is usually called by something else.
Nothing to see here. Move along.
"""
function printhist(tablein, headers; TeX = false, figures = 5)
    ntab = length(tablein[1, :])
    bighead = Array{String,2}(undef, 1, ntab + 1)
    bighead[1] = "n"
    for ih = 2:ntab+1
        bighead[ih] = headers[ih-1]
    end
    if ntab > 5
        error("Too many columns for the table. Use fewer.")
    end
    fmtout = buildformat(ntab, TeX, figures)
    tabfmt = fmtout.format
    nanspace = fmtout.nanspace
    headerfmt = fmtout.headerfmt
    printf(fmt::String, args...) = @eval @printf($fmt, $(args...))
    sprintf(fmt::String, args...) = @eval @sprintf($fmt, $(args...))
    itmax = length(tablein[:, 1])
    itc = 0:itmax-1
    if TeX
        @printf("\\begin{tabular}{")
        for i = 1:ntab+1
            @printf("l")
        end
        @printf("} \n")
        printf(headerfmt, bighead...)
    else
        printf(headerfmt, bighead...)
    end
    for it = 1:itmax
        st = sprintf(tabfmt, it - 1, tablein[it, :]...)
        snan = findfirst(isequal('N'), st)
        lt = length(st)
        while typeof(snan) != Nothing
            st = string(st[1:snan-1], nanspace, st[snan+3:lt])
            snan = findfirst(isequal('N'), st)
        end
        printf("%s", st)
    end
    if TeX
        @printf("\\hline \n")
        @printf("\\end{tabular} \n")
    end
end

function buildformat(ntab, TeX, figures)
    format = "%2d "
    nanspace = "   "
    fwide = figures + 7
    fm1 = string(fwide)
    fm2 = string(figures)
    basefmt = string(" %", fm1, ".", fm2, "e ")
    hfmt = string(" %", fm1, "s")
    headerfmt = " %s"
    if TeX
        for i = 1:ntab
            format = string(format, "&", basefmt)
            headerfmt = string(headerfmt, "&", hfmt)
        end
        format = string(format, " \\\\ \n")
        headerfmt = string(headerfmt, " \\\\ \\hline \n")
    else
        for i = 1:ntab
            format = string(format, basefmt)
            headerfmt = string(headerfmt, hfmt)
        end
        format = string(format, " \n")
        headerfmt = string(headerfmt, " \n \n")
    end
    return (format = format, nanspace = nanspace, headerfmt = headerfmt)
end


text:
module operators_lazytensor

export LazyTensor

import Base: ==, *, /, +, -
import ..operators

using ..sortedindices, ..bases, ..states, ..operators
using ..operators_dense, ..operators_sparse

"""
    LazyTensor(b1[, b2], indices, operators[, factor=1])

Lazy implementation of a tensor product of operators.

The suboperators are stored in the `operators` field. The `indices` field
specifies in which subsystem the corresponding operator lives. Additionally,
a complex factor is stored in the `factor` field which allows for fast
multiplication with numbers.
"""
mutable struct LazyTensor <: Operator
    basis_l::CompositeBasis
    basis_r::CompositeBasis
    factor::Complex128
    indices::Vector{Int}
    operators::Vector{Operator}

    function LazyTensor(op::LazyTensor, factor::Number)
        new(op.basis_l, op.basis_r, factor, op.indices, op.operators)
    end

    function LazyTensor(basis_l::Basis, basis_r::Basis,
                        indices::Vector{Int}, ops::Vector,
                        factor::Number=1)
        if typeof(basis_l) != CompositeBasis
            basis_l = CompositeBasis(basis_l.shape, Basis[basis_l])
        end
        if typeof(basis_r) != CompositeBasis
            basis_r = CompositeBasis(basis_r.shape, Basis[basis_r])
        end
        N = length(basis_l.bases)
        @assert N==length(basis_r.bases)
        sortedindices.check_indices(N, indices)
        @assert length(indices) == length(ops)
        for n=1:length(indices)
            @assert isa(ops[n], Operator)
            @assert ops[n].basis_l == basis_l.bases[indices[n]]
            @assert ops[n].basis_r == basis_r.bases[indices[n]]
        end
        if !issorted(indices)
            perm = sortperm(indices)
            indices = indices[perm]
            ops = ops[perm]
        end
        new(basis_l, basis_r, complex(factor), indices, ops)
    end
end

LazyTensor(basis::Basis, indices::Vector{Int}, ops::Vector, factor::Number=1) = LazyTensor(basis, basis, indices, ops, factor)
LazyTensor(basis_l::Basis, basis_r::Basis, index::Int, operator::Operator, factor::Number=1) = LazyTensor(basis_l, basis_r, [index], Operator[operator], factor)
LazyTensor(basis::Basis, index::Int, operators::Operator, factor::Number=1.) = LazyTensor(basis, basis, index, operators, factor)

Base.copy(x::LazyTensor) = LazyTensor(x.basis_l, x.basis_r, copy(x.indices), [copy(op) for op in x.operators], x.factor)

"""
    suboperator(op::LazyTensor, index)

Return the suboperator corresponding to the subsystem specified by `index`. Fails
if there is no corresponding operator (i.e. it would be an identity operater).
"""
suboperator(op::LazyTensor, index::Int) = op.operators[findfirst(op.indices, index)]
"""
    suboperators(op::LazyTensor, index)

Return the suboperators corresponding to the subsystems specified by `indices`. Fails
if there is no corresponding operator (i.e. it would be an identity operater).
"""
suboperators(op::LazyTensor, indices::Vector{Int}) = op.operators[[findfirst(op.indices, i) for i in indices]]

Base.full(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op.indices, DenseOperator[full(x) for x in op.operators])
Base.sparse(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op.indices, SparseOperator[sparse(x) for x in op.operators])

==(x::LazyTensor, y::LazyTensor) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.operators==y.operators && x.factor==y.factor

# Arithmetic operations
-(a::LazyTensor) = LazyTensor(a, -a.factor)

function *(a::LazyTensor, b::LazyTensor)
    check_multiplicable(a, b)
    indices = sortedindices.union(a.indices, b.indices)
    ops = Vector{Operator}(length(indices))
    for n in 1:length(indices)
        i = indices[n]
        in_a = i in a.indices
        in_b = i in b.indices
        if in_a && in_b
            ops[n] = suboperator(a, i)*suboperator(b, i)
        elseif in_a
            a_i = suboperator(a, i)
            ops[n] = a_i*identityoperator(typeof(a_i), b.basis_l.bases[i], b.basis_r.bases[i])
        elseif in_b
            b_i = suboperator(b, i)
            ops[n] = identityoperator(typeof(b_i), a.basis_l.bases[i], a.basis_r.bases[i])*b_i
        end
    end
    return LazyTensor(a.basis_l, b.basis_r, indices, ops, a.factor*b.factor)
end
*(a::LazyTensor, b::Number) = LazyTensor(a, a.factor*b)
*(a::Number, b::LazyTensor) = LazyTensor(b, a*b.factor)
function *(a::LazyTensor, b::DenseOperator)
    check_multiplicable(a, b)
    result = DenseOperator(a.basis_l, b.basis_r)
    operators.gemm!(complex(1.), a, b, complex(1.), result)
    result
end
function *(a::DenseOperator, b::LazyTensor)
    check_multiplicable(a, b)
    result = DenseOperator(a.basis_l, b.basis_r)
    operators.gemm!(complex(1.), a, b, complex(1.), result)
    result
end

/(a::LazyTensor, b::Number) = LazyTensor(a, a.factor/b)

operators.dagger(op::LazyTensor) = LazyTensor(op.basis_r, op.basis_l, op.indices, Operator[dagger(x) for x in op.operators], conj(op.factor))

operators.tensor(a::LazyTensor, b::LazyTensor) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [a.indices; b.indices+length(a.basis_l.bases)], Operator[a.operators; b.operators], a.factor*b.factor)

function operators.trace(op::LazyTensor)
    b = basis(op)
    result = op.factor
    for i in 1:length(b.bases)
        if i in op.indices
            result *= trace(suboperator(op, i))
        else
            result *= length(b.bases[i])
        end
    end
    result
end

function operators.ptrace(op::LazyTensor, indices::Vector{Int})
    operators.check_ptrace_arguments(op, indices)
    N = length(op.basis_l.shape)
    rank = N - length(indices)
    factor = op.factor
    for i in indices
        if i in op.indices
            factor *= trace(suboperator(op, i))
        else
            factor *= length(op.basis_l.bases[i])
        end
    end
    remaining_indices = sortedindices.remove(op.indices, indices)
    if rank==1 && length(remaining_indices)==1
        return factor * suboperator(op, remaining_indices[1])
    end
    b_l = ptrace(op.basis_l, indices)
    b_r = ptrace(op.basis_r, indices)
    if rank==1
        return factor * identityoperator(b_l, b_r)
    end
    ops = Vector{Operator}(length(remaining_indices))
    for i in 1:length(ops)
        ops[i] = suboperator(op, remaining_indices[i])
    end
    LazyTensor(b_l, b_r, sortedindices.shiftremove(op.indices, indices), ops, factor)
end

operators.normalize!(op::LazyTensor) = (op.factor /= trace(op))

function operators.permutesystems(op::LazyTensor, perm::Vector{Int})
    b_l = permutesystems(op.basis_l, perm)
    b_r = permutesystems(op.basis_r, perm)
    indices = [findfirst(perm, i) for i in op.indices]
    perm_ = sortperm(indices)
    LazyTensor(b_l, b_r, indices[perm_], op.operators[perm_], op.factor)
end

operators.identityoperator(::Type{LazyTensor}, b1::Basis, b2::Basis) = LazyTensor(b1, b2, Int[], Operator[])

# Recursively calculate result_{IK} = \\sum_J op_{IJ} h_{JK}
function _gemm_recursive_dense_lazy(i_k::Int, N_k::Int, K::Int, J::Int, val::Complex128,
                        shape::Vector{Int}, strides_k::Vector{Int}, strides_j::Vector{Int},
                        indices::Vector{Int}, h::LazyTensor,
                        op::Matrix{Complex128}, result::Matrix{Complex128})
    if i_k > N_k
        for I=1:size(op, 1)
            result[I, K] += val*op[I, J]
        end
        return nothing
    end
    if i_k in indices
        h_i = operators_lazytensor.suboperator(h, i_k)
        if isa(h_i, SparseOperator)
            h_i_data = h_i.data::SparseMatrixCSC{Complex128,Int}
            @inbounds for k=1:h_i_data.n
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for jptr=h_i_data.colptr[k]:h_i_data.colptr[k+1]-1
                    j = h_i_data.rowval[jptr]
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data.nzval[jptr]
                    _gemm_recursive_dense_lazy(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        elseif isa(h_i, DenseOperator)
            h_i_data = h_i.data::Matrix{Complex128}
            Nk = size(h_i_data, 2)
            Nj = size(h_i_data, 1)
            @inbounds for k=1:Nk
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for j=1:Nj
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data[j,k]
                    _gemm_recursive_dense_lazy(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        else
            throw(ArgumentError("gemm! of LazyTensor is not implemented for $(typeof(h_i))"))
        end
    else
        @inbounds for k=1:shape[i_k]
            K_ = K + strides_k[i_k]*(k-1)
            J_ = J + strides_j[i_k]*(k-1)
            _gemm_recursive_dense_lazy(i_k + 1, N_k, K_, J_, val, shape, strides_k, strides_j, indices, h, op, result)
        end
    end
end

# Recursively calculate result_{JI} = \\sum_K h_{JK} op_{KI}
function _gemm_recursive_lazy_dense(i_k::Int, N_k::Int, K::Int, J::Int, val::Complex128,
                        shape::Vector{Int}, strides_k::Vector{Int}, strides_j::Vector{Int},
                        indices::Vector{Int}, h::LazyTensor,
                        op::Matrix{Complex128}, result::Matrix{Complex128})
    if i_k > N_k
        for I=1:size(op, 2)
            result[J, I] += val*op[K, I]
        end
        return nothing
    end
    if i_k in indices
        h_i = suboperator(h, i_k)
        if isa(h_i, SparseOperator)
            h_i_data = h_i.data::SparseMatrixCSC{Complex128,Int}
            @inbounds for k=1:h_i_data.n
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for jptr=h_i_data.colptr[k]:h_i_data.colptr[k+1]-1
                    j = h_i_data.rowval[jptr]
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data.nzval[jptr]
                    _gemm_recursive_lazy_dense(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        elseif isa(h_i, DenseOperator)
            h_i_data = h_i.data::Matrix{Complex128}
            Nk = size(h_i_data, 2)
            Nj = size(h_i_data, 1)
            @inbounds for k=1:Nk
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for j=1:Nj
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data[j,k]
                    _gemm_recursive_lazy_dense(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        else
            throw(ArgumentError("gemm! of LazyTensor is not implemented for $(typeof(h_i))"))
        end
    else
        @inbounds for k=1:shape[i_k]
            K_ = K + strides_k[i_k]*(k-1)
            J_ = J + strides_j[i_k]*(k-1)
            _gemm_recursive_lazy_dense(i_k + 1, N_k, K_, J_, val, shape, strides_k, strides_j, indices, h, op, result)
        end
    end
end

function gemm(alpha::Complex128, op::Matrix{Complex128}, h::LazyTensor, beta::Complex128, result::Matrix{Complex128})
    if beta == Complex128(0.)
        fill!(result, beta)
    elseif beta != Complex128(1.)
        scale!(beta, result)
    end
    N_k = length(h.basis_r.bases)
    shape = [min(h.basis_l.shape[i], h.basis_r.shape[i]) for i=1:length(h.basis_l.shape)]
    strides_j = operators_dense._strides(h.basis_l.shape)
    strides_k = operators_dense._strides(h.basis_r.shape)
    _gemm_recursive_dense_lazy(1, N_k, 1, 1, alpha*h.factor, shape, strides_k, strides_j, h.indices, h, op, result)
end

function gemm(alpha::Complex128, h::LazyTensor, op::Matrix{Complex128}, beta::Complex128, result::Matrix{Complex128})
    if beta == Complex128(0.)
        fill!(result, beta)
    elseif beta != Complex128(1.)
        scale!(beta, result)
    end
    N_k = length(h.basis_l.bases)
    shape = [min(h.basis_l.shape[i], h.basis_r.shape[i]) for i=1:length(h.basis_l.shape)]
    strides_j = operators_dense._strides(h.basis_l.shape)
    strides_k = operators_dense._strides(h.basis_r.shape)
    _gemm_recursive_lazy_dense(1, N_k, 1, 1, alpha*h.factor, shape, strides_k, strides_j, h.indices, h, op, result)
end

operators.gemm!(alpha, h::LazyTensor, op::DenseOperator, beta, result::DenseOperator) = gemm(convert(Complex128, alpha), h, op.data, convert(Complex128, beta), result.data)
operators.gemm!(alpha, op::DenseOperator, h::LazyTensor, beta, result::DenseOperator) = gemm(convert(Complex128, alpha), op.data, h, convert(Complex128, beta), result.data)

function operators.gemv!(alpha::Complex128, a::LazyTensor, b::Ket, beta::Complex128, result::Ket)
    b_data = reshape(b.data, length(b.data), 1)
    result_data = reshape(result.data, length(result.data), 1)
    gemm(alpha, a, b_data, beta, result_data)
end

function operators.gemv!(alpha::Complex128, a::Bra, b::LazyTensor, beta::Complex128, result::Bra)
    a_data = reshape(a.data, 1, length(a.data))
    result_data = reshape(result.data, 1, length(result.data))
    gemm(alpha, a_data, b, beta, result_data)
end

operators.gemv!(alpha, a::LazyTensor, b::Ket, beta, result::Ket) = operators.gemv!(convert(Complex128, alpha), a, b, convert(Complex128, beta), result)
operators.gemv!(alpha, a::Bra, b::LazyTensor, beta, result::Bra) = operators.gemv!(convert(Complex128, alpha), a, b, convert(Complex128, beta), result)

end # module


text:
function initialize_pss!(
    device_states,
    device::PSY.DynamicGenerator{M, S, A, TG, PSY.PSSFixed},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov}

end

#Currently not working properly.
#=
function initialize_pss!(
    device_states,
    device::PSY.DynamicGenerator{M, S, A, TG, PSY.PSSSimple},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov}

    ω0 = PSY.get_ext(device)[CONTROL_REFS][ω_ref_index]
    static_gen = PSY.get_static_injector(device)
    P0 = PSY.get_active_power(static_gen)
    τe = get_inner_vars(device)[τe_var]

    #Get parameters
    pss = PSY.get_pss(device)
    K_ω = PSY.get_K_ω(pss)
    K_p = PSY.get_K_p(pss)

    #It assumes that the system is balanced so:
    get_inner_vars(device)[V_pss_var] = K_ω * (ω0 - ω0) + K_p * (ω0 * τe - P0)

end
=#


text:
# Advent of Code 2021, Day 15
#
# Objective is to start at the top right, and end up at
# the bottom right, adding up any cells you enter, so as to
# minimize total cost.
#
# First attempt depth-first search with maximum look-ahead,
# but did not succeed on full data set. New solution uses
# Dijkstra's algorithm, adapted from:
# https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
# Note that the numbered steps in the comments are copied
# from this Wikipedia article.
#
# AK, 15/12/2021

# Name of input file
f = "sample.txt"
f = "input.txt"

# Read data into a matrix
d = split(strip(read(f, String)), "\n")
l = map(r -> [x - '0' for x in r], d)
m = Matrix(transpose(reshape(vcat(l...), (:, length(l)))))
#display(m)

# Use this value to represent extremely high cost
PENALTY = 99999999

# Find the lowest total cost path from 1,1 to m,n of a matrix
function solve(M::Matrix{Int64})

    # 1. Mark all nodes as unvisited
    nr, nc = size(M)
    visited = zeros(nr, nc)

    # 2. Assign to every node a tentative distance value: set it to zero for our
    # initial node and to infinity for all other nodes. The tentative distance of a
    # node v is the length of the shortest path discovered so far between the node
    # v and the starting node. Since initially no path is known to any other vertex
    # than the source itself (which is a path of length zero), all other tentative
    # distances are initially set to infinity. Set the initial node as current.[15]
    distances = fill(PENALTY, nr, nc)
    distances[1,1] = 0

    # Start algoritm iterations, will stop when reach bottom right
    iter = 0
    r, c = 1, 1     # Initial node: top left
    while ! (r == nr && c == nc)

        # 3. For the current node, consider all of its unvisited neighbors and
        # calculate their tentative distances through the current node. Compare the
        # newly calculated tentative distance to the current assigned value and assign
        # the smaller one. For example, if the current node A is marked with a distance
        # of 6, and the edge connecting it with a neighbor B has length 2, then the
        # distance to B through A will be 6 + 2 = 8. If B was previously marked with a
        # distance greater than 8 then change it to 8. Otherwise, the current value
        # will be kept.
        for p in [[r, c-1], [r, c+1], [r-1,c], [r+1,c]] # left, right, up, down
            y, x = p
            if y >= 1 && y <= nr && x >= 1 && x <= nc && visited[y, x] == 0
                c1 = M[y, x] + distances[r, c]
                if c1 < distances[y, x]
                    distances[y, x] = c1
                end
            end
        end

        # 4. When we are done considering all of the unvisited neighbors of the current
        # node, mark the current node as visited and remove it from the unvisited set.
        # A visited node will never be checked again.
        visited[r, c] = 1

        # 5. If the destination node has been marked visited (when planning a route
        # between two specific nodes) or if the smallest tentative distance among the
        # nodes in the unvisited set is infinity (when planning a complete traversal;
        # occurs when there is no connection between the initial node and remaining
        # unvisited nodes), then stop. The algorithm has finished.
        #
        # 6. Otherwise, select the unvisited node that is marked with the smallest
        # tentative distance, set it as the new current node, and go back to step 3.
        #
        # When planning a route, it is actually not necessary to wait until the
        # destination node is "visited" as above: the algorithm can stop once the
        # destination node has the smallest tentative distance among all "unvisited"
        # nodes (and thus could be selected as the next "current").
        #
        # In this case, get the coordinates of the minimum value in the matrix, first
        # changing visited cells to high-cost penalty so they are not considered.
        nextCosts = distances .+ visited * PENALTY
        r, c = argmin(nextCosts).I   # coords of min cost (next iter starts there)

        # Report progress
        iter += 1
        if iter % 1000 == 0
            print("\rIteration $iter: now at $r,$c   ")
        end

    end

    # Return solution, cost set in the distances matrix for bottom right point
    return distances[nr, nc]

end

# For part 2, make a 5x5 copy of the matrix, each row/col being one higher
# than the previous, with 9s wrapped back to 1
function enlarge(x, n)

    # To grow 5x, add 4 copies either direction
    n -= 1

    # Grow horizontally
    m1 = copy(x)
    for i in 1:n
        m1 = hcat(m1, x .+ i)
    end

    # Grow vertically
    m2 = copy(m1)
    for i in 1:n
        m2 = vcat(m2, m1 .+ i)
    end

    # Adjust cells that are bigger than 9
    nr, nc = size(m2)
    for r in 1:nr
        for c in 1:nc
            while m2[r,c] > 9
                m2[r,c] -= 9
            end
        end
    end

    # Return result
    return m2
end

# Part 1 solution
s1 = solve(m)
println("\nPart 1: ", s1)

# Part 2, same but on enlarged matrix
mbig = enlarge(m, 5)
s2 = solve(mbig)
println("\nPart 2: ", s2)


text:
using QuasiGeostrophy

abstract type SpatialDiscretization <: AbstractExpression end
abstract type TimeDiscretization <: AbstractExpression end
struct Fourier <: SpatialDiscretization end
struct Time <: TimeDiscretization end

struct AB1 <: TimeDiscretization end
struct BDF1 <: TimeDiscretization end

struct Explicit{T,S} <: TimeDiscretization
    operand::T
    meta_data::S
end

function extract(a, b)
    if typeof(a) <: b
        return a
    else
        return Field(0)
    end
    return nothing
end

##
struct Implicit{T,S} <: TimeDiscretization
    operand::T
    meta_data::S
end

struct AbstractEquation{T,S} <: AbstractExpression
    lhs::T
    rhs::S
end
struct Wrapper{T, S} <: AbstractExpression
    data::T
    meta_data::S
end
# Struct for MetaData
struct MetaData{T}
    io_name::T
end
struct GradientMetaData{T}
    direction::T
end
function Base.show(io::IO, field::Wrapper{T, S}) where {T, S <: MetaData}
    color = 226
    printstyled(io, field.meta_data.io_name, color = color)
end

u = Wrapper("u", MetaData("u"))
σ = Wrapper("σ", MetaData("σ"))
ν = Wrapper(0.001, MetaData("ν"))
∂x(u) = Gradient(u, GradientMetaData("x"))
∂t(u) = Gradient(u, GradientMetaData("t"))
tmp = ∂x(u)
##
function Base.show(io::IO, ∇::Gradient)
    color = 220
    printstyled(io, "∇(", color = color)
    print(∇.operand)
    printstyled(io, ")", color = color)
end

function Base.show(io::IO, ∇::Gradient)
    color = 220
    printstyled(io, "∇(", color = color)
    print(∇.operand)
    printstyled(io, ")", color = color)
end

function Base.show(io::IO, ∇::Gradient{T,S}) where {T, S <: GradientMetaData}
    color = 220
    tmp = ∇.metadata.direction
    printstyled(io, "∂", tmp, "(", color = color)
    print(∇.operand)
    printstyled(io, ")", color = color)
end

function Base.show(io::IO, ϕ::Field)
    color = 220
    print(ϕ.data)
end
##
import Base: ==

function ==(a::AbstractExpression, b::AbstractExpression)
    return AbstractEquation(a,b)
end

function Base.show(io::IO, eq::AbstractEquation)
    color = 118
    print(io, eq.lhs)
    printstyled(io, "=", color = color)
    print(io, eq.rhs)
end

function Base.show(io::IO, eq::Implicit)
    print(io, eq.operand)
end
function Base.show(io::IO, eq::Explicit)
    print(io, eq.operand)
end

function Base.show(io::IO, system::Array{T}) where {T <: AbstractEquation}
    for equation in system
        println(io, equation)
    end
end

∂t(u) == -∂x(u*u) + ν * ∂x(∂x(u))

pde_equation = [
    Implicit(∂t(u), BDF1()) == Explicit(-∂x(u*u), AB1()) + Implicit(ν * ∂x(∂x(u)) , BDF1())
]

## Implicit / Explicit grabbing rules.
# Assumes seperation by addition should be okay since quadrature is a
# linear operation
function get_implicit(a::Implicit)
    return a
end
function get_implicit(a::Explicit)
    return Field(0) # probably need a rule like 0 + field = field
end
function get_implicit(a::Add)
    t1 = get_implicit(a.term1)
    t2 = get_implicit(a.term2)
    return t1 + t2
end
function get_implicit(a::Negative)
    t1 = get_implicit(a.term)
    return t1
end

function get_explicit(a::Explicit)
    return a
end

function get_explicit(a::Implicit)
    return Field(0) # probably need a rule like 0 + field = field
end

function get_explicit(a::Add)
    t1 = get_explicit(a.term1)
    t2 = get_explicit(a.term2)
    return t1 + t2
end
function get_explicit(a::Negative)
    t1 = get_explicit(a.term)
    return t1
end

function get_implicit(eq::AbstractEquation)
    # place all implicit terms on left-hand side
    implicit = get_implicit(eq.lhs) - get_implicit(eq.rhs)
    return implicit
end
function get_explicit(eq::AbstractEquation)
    # place all explicit terms on left-hand side
    explicit = get_explicit(eq.rhs) - get_explicit(eq.lhs)
    return explicit
end
function get_imex(eq::AbstractEquation)
    implicit = get_implicit(eq)
    explicit = get_explicit(eq)
    return implicit == explicit
end
function get_imex(system::Array{S}) where {S <: AbstractEquation}
    return [get_imex(equation) for equation in system]
end
##
tmp_expr = Explicit(-∂x(u*u), AB1()) + Implicit(ν * ∂x(∂x(u)) , BDF1()) + Explicit(-∂x(u*u), AB1()) + Explicit(-∂x(u*u), AB1())
tmp_implicit = get_implicit(tmp_expr)
tmp_implicit = get_explicit(tmp_expr)
new_eq  = get_imex(Implicit(∂t(u), BDF1()) == Explicit(-∂x(u*u), AB1()) + Implicit(ν * ∂x(∂x(u)) , BDF1()))
##
pde_system = [
    Implicit(σ, BDF1()) == Implicit(∂x(u), BDF1()),
    Implicit(∂t(u), BDF1()) == Explicit(- ∂x( u * u), AB1())+ Implicit(ν *  ∂x(σ), BDF1())
]
pde_system[1]
pde_system[2]
new_system = get_imex(pde_system)
eq1 = pde_system[1]
eq2 = pde_system[2]
println("The system of equations are")
for equation in pde_system
    println(equation)
end
println("The new system with all implicit terms
         on the left and explicit on the right is")
for equation in new_system
    println(equation)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
#========== backward gradient using ForwardDiff ==========#

function insert_forward_gradient(axislist, store)
    dZ = Symbol(DEL, ZED)
    ∇act! = Symbol(:∇, ACT!)
    gradarrays = map(A -> Symbol(DEL, A), store.arrays)
    # gradscalars = map(A -> Symbol(DEL, A), store.scalars)
    nonshared = setdiff(vcat(store.leftind, store.redind), store.sharedind)

    epsilondict = Dict{Symbol,Expr}()

    epsilonright = MacroTools_postwalk(epsilonwalk(epsilondict), store.right)
    # epsilonright = MacroTools_postwalk(epsilonwalk(epsilondict, store.scalars), store.right)

    defineepsilons, readepsilons = [], []
    for (d, (Aepsilon, Aex)) in enumerate(epsilondict)
        basis = [i==d ? :(one($TYP)) : :(zero($TYP)) for i in 1:length(epsilondict)]
        push!(defineepsilons, :($Aepsilon = ForwardDiff.Dual(zero($TYP), ($(basis...),))))
        push!(readepsilons, :($Aex = $Aex + ForwardDiff.partials($ZED, $d) * $dZ[$(store.leftraw...)]))
    end

    ex_iter = :($ZED = $(epsilonright); $(readepsilons...))

    make_many_actors(∇act!,
        vcat(gradarrays, :($dZ::AbstractArray{$TYP}), ZED, store.arrays, store.scalars, axislist),
        # vcat(gradarrays, gradscalars, :($dZ::AbstractArray{$TYP}), store.arrays, store.scalars, axislist),
        :(($(defineepsilons...);)), store.sharedind, nothing, nonshared, ex_iter, nothing, store, " (gradient using ForwardDiff)")

    if isdefined(store.mod, :Zygote)
        ex_iter2 = fillarrayreplace(ex_iter, dZ)
        ex_value = :($(Symbol(dZ, :_value)) = $dZ.value)

        make_many_actors(∇act!,
            vcat(gradarrays, :($dZ::Zygote.Fill{$TYP}), ZED, store.arrays, store.scalars, axislist),
            :(($(defineepsilons...); $ex_value)), store.sharedind, nothing, nonshared, ex_iter2, nothing, store, " (method for FillArrays)")

        # push!(store.outeval, quote
        #     Tullio.promote_storage(T::Type, ::Type{<:Zygote.Fill}) = T
        #     Tullio.promote_storage(::Type{<:Zygote.Fill}, T::Type) = T
        # end)
    end

end

epsilonwalk(dict) = ex -> begin
# epsilonwalk(dict, scalars) = ex -> begin
#         ex isa Symbol && ex in scalars && return scalarplusepsilon(ex, dict)
        @capture_(ex, A_[inds__]) || return ex
        return arrayplusepsilon(A, inds, dict)
    end

arrayplusepsilon(A::Symbol, inds, dict) = begin # the same array may occur twice!
    Aepsilon = Symbol(EPS, A)
    while haskey(dict, Aepsilon)
        Aepsilon = Symbol(Aepsilon, "′")
    end
    dict[Aepsilon] = :( $(Symbol(DEL, A))[$(inds...)] )
    :(( $A[$(inds...)] + $Aepsilon ))
end
arrayplusepsilon(A, inds, dict) = begin
    @debug string("expression ", A, " is why you can't use ForwardDiff here")
    :🐳
end

# scalarplusepsilon(A::Symbol, dict) = begin
#     Aepsilon = Symbol(EPS, A)
#     dict[Aepsilon] = Symbol(DEL, A)
#     :(( $A + $Aepsilon ))
# end

#========== the end ==========#


text:
export reidentify
export relevance_mask

"""
    reidentify(z::AbstractMatrix, zmask::AbstractMatrix,
               u::AbstractMatrix, decoder::ACD, opt; N=200)

Optimize MSE(decoder(z), u) where only the values of z that are specified by
zmask are allowed to change. z and u are assumed to have samples in the second
dimension.

# Arguments
* `z`: latent variables typically sampled from the encoder of an AbstractGM
* `zmask`: specifies that values of z that are allowed to change during optimisation
* `u`: labels for optimisation
* `decoder`: decodes z into u-space
* `opt`: Flux optimiser
* `N`: optimisation steps
"""
function reidentify(z::AbstractMatrix, zmask::AbstractMatrix, u::AbstractMatrix,
                    decoder::ACD, opt; N=200)
    zr = copy(z)
    ps = Flux.Params([zr])
    loss() = -sum(loglikelihood(decoder, u, zr))

    for ii in 1:N
        try
            d = Dict()
            gs = Flux.gradient(loss, ps)
            for p in ps
                g = gs[p]
                g[zmask] .= 0
                d[p] = g
            end
            gs = Zygote.Grads(d)
            Flux.Optimise.update!(opt, ps, gs)
            @info "ii=$ii N=$N loss=$(loss())"
        catch ex
            rethrow(ex)
        end
    end
    return zr
end

"""
    reidentify(zmask::AbstractVector, u::AbstractVector, zbatchsize::Int,
               model::AbstractGM, opt; N=100)

Optimize MSE(decoder(z), u) where only values of z that are specified by zmask
are allowed to change. z is sampled `zbatchsize` times via
`rand(model.encoder,u)`. After `N` steps with the given optimizer only the best
of the sampled `z`s is returned.
"""
function reidentify(zmask::AbstractVector, u::AbstractVector, zbatchsize::Int,
                    model::AbstractGM, opt; N=100)
    U = repeat(u, inner=(1, zbatchsize))
    Z = rand(model.encoder, U)
    X = mean(model.decoder, Z)

    zlen = size(Z,1)
    Zmask = repeat(reshape(.!zmask, zlen, 1), inner=(1, zbatchsize))
    Zr = reidentify(Z, Zmask, U, model.decoder, opt, N=N)

    # get best result
    llh = dropdims(loglikelihood(model.decoder, U, Zr), dims=1)
    idx = argmax(llh)
    z   = Z[:,idx]
    x   = X[:,idx]
    zr  = Zr[:,idx]
    xr  = mean(model.decoder, zr)

    (z, x, zr, xr)
end

function relevance_mask(Z::AbstractMatrix; σmax=0.05, μmax=0.05)
    μ = dropdims(mean(Z, dims=2), dims=2)
    σ = dropdims(std(Z, dims=2), dims=2)
    s = σ .> σmax
    m = abs.(μ) .> μmax
    m .| s
end


text:
using BasicTextRender, FileIO, ImageCore

@time img = rendertext("Hello world", height=50)
save(joinpath(@__DIR__,"img.png"), img)

img = rendertext("Hello world", height=20, color=RGBA(1.0,0.0,0.0,1.0), backgroundColor=RGBA(1.0,1.0,0.0,1.0))
save(joinpath(@__DIR__,"img2.png"), img)

img = rand(RGBA{Float64},80,700)
overlaytext!(img, "BasicTextRender.jl", 60, (10,40))
save(joinpath(@__DIR__,"logo.png"), img)

using BenchmarkTools
@btime rendertext("Hello world", height=30)


text:
using KernelAbstractions
using Test
using CUDAapi
if has_cuda_gpu()
    using CuArrays
    CuArrays.allowscalar(false)
end

@kernel function private(A)
    @uniform N = prod(groupsize())
    I = @index(Global, Linear)
    i = @index(Local, Linear)
    priv = @private Int (1,)
    @inbounds begin
        priv[1] = N - i + 1
        @synchronize
        A[I] = priv[1]
    end
end

# This is horrible don't write code like this
@kernel function forloop(A, ::Val{N}) where N
    I = @index(Global, Linear)
    i = @index(Local, Linear)
    priv = @private Int (N,)
    for j in 1:N
        priv[j] = A[I, j]
    end
    A[I, 1] = 0
    @synchronize
    for j in 1:N
        k = mod1(j + i - 1, N)
        A[k, 1] += priv[j]
        @synchronize
    end
end

function harness(backend, ArrayT)
    A = ArrayT{Int}(undef, 64)
    wait(private(backend, 16)(A, ndrange=size(A)))
    @test all(A[1:16] .== 16:-1:1)
    @test all(A[17:32] .== 16:-1:1)
    @test all(A[33:48] .== 16:-1:1)
    @test all(A[49:64] .== 16:-1:1)

    A = ArrayT{Int}(undef, 64, 64)
    A .= 1
    wait(forloop(backend)(A, Val(size(A, 2)), ndrange=size(A,1), workgroupsize=size(A,1)))
    @test all(A[:, 1] .== 64)
    @test all(A[:, 2:end] .== 1)
end

@testset "kernels" begin
    harness(CPU(), Array)
    if has_cuda_gpu()
        harness(CUDA(), CuArray)
    end
end


text:
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

import Hafta
import HaftaViz

@test HaftaViz.Hafta === Hafta


text:
struct VideoStream
    io
    process
    screen
    path::String
end

"""
    VideoStream(scene::Scene, dir = mktempdir(), name = "video")

returns a stream and a buffer that you can use to not allocate for new frames.
Use `add_frame!(stream, window, buffer)` to add new video frames to the stream.
Use `finish(stream)` to save the video to 'dir/name.mkv'. You can also call
`finish(stream, "mkv")`, `finish(stream, "mp4")`, `finish(stream, "gif")` or `finish(stream, "webm")` to convert the stream to those formats.
"""
function VideoStream(scene::Scene)
    if !has_ffmpeg[]
        error("You can't create a video stream without ffmpeg installed.
         Please install ffmpeg, e.g. via https://ffmpeg.org/download.html.
         When you download the binaries, please make sure that you add the path to your PATH
         environment variable.
         On unix you can install ffmpeg with `sudo apt-get install ffmpeg`.
        ")
    end
    #codec = `-codec:v libvpx -quality good -cpu-used 0 -b:v 500k -qmin 10 -qmax 42 -maxrate 500k -bufsize 1000k -threads 8`
    dir = mktempdir()
    path = joinpath(dir, "$(gensym(:video)).mkv")

    display(scene)
    AbstractPlotting.update!(scene)
    _xdim, _ydim = widths(pixelarea(scene)[])
    xdim = _xdim % 2 == 0 ? _xdim : _xdim + 1
    ydim = _ydim % 2 == 0 ? _ydim : _ydim + 1
    process = open(`ffmpeg -loglevel quiet -f rawvideo -pixel_format rgb24 -r 24 -s:v $(xdim)x$(ydim) -i pipe:0 -vf vflip -y $path`, "w")
    VideoStream(process.in, process, GLMakie.global_gl_screen(), abspath(path))
end

"""
Adds a video frame to the VideoStream
"""
function recordframe!(io::VideoStream)
    #codec = `-codec:v libvpx -quality good -cpu-used 0 -b:v 500k -qmin 10 -qmax 42 -maxrate 500k -bufsize 1000k -threads 8`
    frame = GLMakie.colorbuffer(io.screen)
    _xdim, _ydim = size(frame)
    xdim = _xdim % 2 == 0 ? _xdim : _xdim + 1
    ydim = _ydim % 2 == 0 ? _ydim : _ydim + 1
    frame_out = fill(RGB{N0f8}(1, 1, 1), ydim, xdim)
    for x in 1:_xdim, y in 1:_ydim
        c = frame[(_xdim + 1) - x, y]
        frame_out[y, x] = RGB{N0f8}(c)
    end
    write(io.io, frame_out)
    return
end

"""
    save(path::String, io::VideoStream)

Flushes the video stream and converts the file to the extension found in `path` which can
be `mkv` is default and doesn't need convert, `gif`, `mp4` and `webm`.
`mp4` is recommended for the internet, since it's the most supported format.
`webm` yields the smallest file size, `mp4` and `mk4` are marginally bigger and `gif`s are up to
6 times bigger with same quality!
"""
function save(path::String, io::VideoStream)
    close(io.process)
    wait(io.process)
    p, typ = splitext(path)
    if typ == ".mkv"
        cp(io.path, out)
    elseif typ == ".mp4"
        run(`ffmpeg -loglevel quiet -i $(io.path) -c:v libx264 -preset slow -crf 24 -pix_fmt yuv420p -c:a libvo_aacenc -b:a 128k -y $path`)
    elseif typ == ".webm"
        run(`ffmpeg -loglevel quiet -i $(io.path) -c:v libvpx-vp9 -threads 16 -b:v 2000k -c:a libvorbis -threads 16 -vf scale=iw:ih -y $path`)
    elseif typ == ".gif"
        filters = "fps=15,scale=iw:ih:flags=lanczos"
        palette_path = dirname(io.path)
        pname = joinpath(palette_path, "palette.bmp")
        isfile(pname) && rm(pname, force = true)
        run(`ffmpeg -loglevel quiet -i $(io.path) -vf "$filters,palettegen" -y $pname`)
        run(`ffmpeg -loglevel quiet -i $(io.path) -i $pname -lavfi "$filters [x]; [x][1:v] paletteuse" -y $path`)
        rm(pname, force = true)
    else
        rm(io.path)
        error("Video type $typ not known")
    end
    rm(io.path)
    return path
end

"""
    record(func, scene, path)
usage:
```example
    record(scene, "test.gif") do io
        for i = 1:100
            scene.plots[:color] = ...# animate scene
            recordframe!(io) # record a new frame
        end
    end
```
"""
function record(func, scene, path)
    io = VideoStream(scene)
    func(io)
    save(path, io)
end

"""
    record(func, scene, path, iter)
usage:
```example
    record(scene, "test.gif", 1:100) do i
        scene.plots[:color] = ...# animate scene
    end
```
"""
function record(func, scene, path, iter)
    io = VideoStream(scene)
    for i in iter
        t1 = time()
        func(i)
        recordframe!(io)
        diff = (1/24) - (time() - t1)
        if diff > 0.0
            sleep(diff)
        else
            yield()
        end
    end
    save(path, io)
end

function show(io::IO, mime::MIME"text/html", vs::VideoStream)
    mktempdir() do dir
        path = finish(vs, joinpath(dir, "video.mp4"))
        print(
            io,
            """<video autoplay controls><source src="data:video/x-m4v;base64,""",
            base64encode(open(read, path)),
            """" type="video/mp4"></video>"""
        )
    end
end


text:
using ConvexTests, ProxSDP

@info "Starting ProxSDP tests"
do_tests("ProxSDP", () -> ProxSDP.Optimizer(tol_primal=1e-6, tol_dual=1e-6); exclude = [r"mip", r"exp"], description="Tests run with `tol_primal=1e-6` and `tol_dual=1e-6`.")


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f7a6d7c3-37b9-437d-8b8e-853665ddbae3
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ 65780f00-ed6b-11ea-1ecf-8b35523a7ac0
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			Pkg.PackageSpec(name="Images", version="0.22.4"),
			Pkg.PackageSpec(name="ImageMagick", version="0.7"),
			Pkg.PackageSpec(name="PlutoUI", version="0.7"),
			Pkg.PackageSpec(name="HypertextLiteral", version="0.5"),
			Pkg.PackageSpec(name="OffsetArrays"),
			])

	using Images
	using PlutoUI
	using HypertextLiteral
	using OffsetArrays
end

# ╔═╡ 83eb9ca0-ed68-11ea-0bc5-99a09c68f867
md"_homework 2, version 1_"

# ╔═╡ ac8ff080-ed61-11ea-3650-d9df06123e1f
md"""

# **Homework 2** - _convolutions_
`18.S191`, Spring 2021

`Due Date`: **Friday Mar 5, 2021 at 11:59pm EST**

This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!
"""

# ╔═╡ 911ccbce-ed68-11ea-3606-0384e7580d7c
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Vidhi Katkoria", kerberos_id = "vidhi")

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

# you might need to wait until all other cells in this notebook have completed running.
# scroll down the page to see what's up

# ╔═╡ 8ef13896-ed68-11ea-160b-3550eeabbd7d
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 5f95e01a-ee0a-11ea-030c-9dba276aba92
md"_Let's create a package environment:_"

# ╔═╡ e08781fa-ed61-11ea-13ae-91a49b5eb74a
md"""

## **Exercise 1** - _Convolutions in 1D_

As we have seen in the lectures, we can produce cool effects using the mathematical technique of **convolutions**. We input one image $M$ and get a new image $M'$ back.

Conceptually we think of $M$ as a matrix. In practice, in Julia it will be a `Matrix` of color objects, and we may need to take that into account. Ideally, however, we should write a **generic** function that will work for any type of data contained in the matrix.

A convolution works on a small **window** of an image, i.e. a region centered around a given point $(i, j)$. We will suppose that the window is a square region with odd side length $2\ell + 1$, running from $-\ell, \ldots, 0, \ldots, \ell$.

The result of the convolution over a given window, centred at the point $(i, j)$ is a *single number*; this number is the value that we will use for $M'_{i, j}$.
(Note that neighbouring windows overlap.)

To get started, in Exercise 1 we'll restrict ourselves to convolutions in 1D.
So a window is just a 1D region from $-\ell$ to $\ell$.

"""

# ╔═╡ 80108d80-ee09-11ea-0368-31546eb0d3cc
md"""
#### Exercise 1.1

Let's create a vector `v` of random numbers of length `n=100`.
"""

# ╔═╡ 7fcd6230-ee09-11ea-314f-a542d00d582e
n = 100

# ╔═╡ 7fdb34dc-ee09-11ea-366b-ffe10d1aa845
v = rand(n)

# ╔═╡ 7fe9153e-ee09-11ea-15b3-6f24fcc20734
md"_Feel free to experiment with different values!_

Let's use the function `colored_line` to view this 1D number array as a 1D image.
"

# ╔═╡ ff70782e-e8d2-4281-9b24-d45c925f55e2
begin
	colored_line(x::Vector) = hcat(Gray.(Float64.(x)))'
	colored_line(x::Any) = nothing
end

# ╔═╡ 01070e28-ee0f-11ea-1928-a7919d452bdd
colored_line(v)

# ╔═╡ 7522f81e-ee1c-11ea-35af-a17eb257ff1a
md"👉 Try changing `n` and `v` around. Notice that you can run the cell `v = rand(n)` again to regenerate new random values."

# ╔═╡ 801d90c0-ee09-11ea-28d6-61b806de26dc
md"""
#### Exercise 1.2
We need to decide how to handle the **boundary conditions**, i.e. what happens if we try to access a position in the vector `v` beyond `1:n`.  The simplest solution is to assume that $v_{i}$ is 0 outside the original vector; however, this may lead to strange boundary effects.

A better solution is to use the *closest* value that is inside the vector. Effectively we are extending the vector and copying the extreme values into the extended positions. (Indeed, this is one way we could implement this; these extra positions are called **ghost cells**.)

👉 Write a function `extend(v, i)` that checks whether the position $i$ is inside `1:n`. If so, return the $(HTML("<br>")) ``i``th component of `v`; otherwise, return the nearest end value.
"""

# ╔═╡ 802bec56-ee09-11ea-043e-51cf1db02a34
function extend(v::AbstractVector, i)
	println(i)
	if i<1
		return v[1]
	end
	if i>length(v)
		return v[length(v)]
	end
	return v[i]
end

# ╔═╡ b7f3994c-ee1b-11ea-211a-d144db8eafc2
md"_Some test cases:_"

# ╔═╡ 3492b164-7065-48e8-978b-6c96b965d376
example_vector = [0.8, 0.2, 0.1, 0.7, 0.6, 0.4]

# ╔═╡ 02123165-2a0a-49a8-b7a9-458955523511
colored_line(example_vector)

# ╔═╡ 806e5766-ee0f-11ea-1efc-d753cd83d086
md"- Extended with 0:"

# ╔═╡ 38da843a-ee0f-11ea-01df-bfa8b1317d36
colored_line([0, 0, example_vector..., 0, 0])

# ╔═╡ 9bde9f92-ee0f-11ea-27f8-ffef5fce2b3c
md"- Extended with your `extend` function:"

# ╔═╡ 431ba330-0f72-416a-92e9-55f51ff3bcd1
md"""
#### Exercise 1.3
👉 Write (or copy) the `mean` function from Homework 1, which takes a vector and returns the mean.

"""

# ╔═╡ 5fdc5d0d-a52c-476e-b3b5-3b6364b706e4
function mean(v)
	return sum(v)/length(v)
end

# ╔═╡ e84c9cc2-e6e1-46f1-bf4e-9605da5e6f4a
md"""

👉 Write a function `box_blur(v, l)` that blurs a vector `v` with a window of length `l` by averaging the elements within a window from $-\ell$ to $\ell$. This is called a **box blur**. Use your function `extend` to handle the boundaries correctly.

Return a vector of the same size as `v`.
"""

# ╔═╡ ca91bf1a-80ed-11eb-3be0-cd60c9adae9d
colored_line(example_vector)

# ╔═╡ 809f5330-ee09-11ea-0e5b-415044b6ac1f
md"""
#### Exercise 1.4
👉 Apply the box blur to your vector `v`. Show the original and the new vector by creating two cells that call `colored_line`. Make the parameter $\ell$ interactive, and call it `l_box` instead of `l` to avoid a naming conflict.
"""

# ╔═╡ e555a7e6-f11a-43ac-8218-6d832f0ce251
@bind l_box Slider(1:20)

# ╔═╡ 302f0842-453f-47bd-a74c-7942d8c96485
colored_line(v)

# ╔═╡ 80ab64f4-ee09-11ea-29b4-498112ed0799
md"""
#### Exercise 1.5
The box blur is a simple example of a **convolution**, i.e. a linear function of a window around each point, given by

$$v'_{i} = \sum_{m}  \, v_{i - m} \, k_{m},$$

where $k$ is a vector called a **kernel**.

Again, we need to take care about what happens if $v_{i -m }$ falls off the end of the vector.

👉 Write a function `convolve(v, k)` that performs this convolution. You need to think of the vector $k$ as being *centred* on the position $i$. So $m$ in the above formula runs between $-\ell$ and $\ell$, where $2\ell + 1$ is the length of the vector $k$.

   You will either need to do the necessary manipulation of indices by hand, or use the `OffsetArrays.jl` package.
"""

# ╔═╡ cf73f9f8-ee12-11ea-39ae-0107e9107ef5
md"_Edit the cell above, or create a new cell with your own test cases!_"

# ╔═╡ fa463b71-5aa4-44a3-a67b-6b0776236243
md"""
#### Exercise 1.6

👉 Define a function `box_blur_kernel(l)` which returns a _kernel_ (i.e. a vector) which, when used as the kernel in `convolve`, will reproduce a box blur of length `l`.
"""

# ╔═╡ 8a7d3cfd-6f19-43f0-ae16-d5a236f148e7
function box_blur_kernel(l)

	kernel_size = 2*l+1
	box_blur_kernel = fill(1/kernel_size, kernel_size)
	return box_blur_kernel
end

# ╔═╡ a34d1ad8-3776-4bc4-93e5-72cfffc54f15
@bind box_kernel_l Slider(1:5)

# ╔═╡ 971a801d-9c46-417a-ad31-1144894fb4e1
box_blur_kernel_test = box_blur_kernel(box_kernel_l)

# ╔═╡ 5f13b1a5-8c7d-47c9-b96a-a09faf38fe5e
md"""
Let's apply your kernel to our test vector `v` (first cell), and compare the result to our previous box blur function (second cell). The two should be identical.
"""

# ╔═╡ 03f91a22-1c3e-4c42-9d78-1ee36851a120
md"""
#### Exercise 1.7
👉 Write a function `gaussian_kernel`.

The definition of a Gaussian in 1D is

$$G(x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( \frac{-x^2}{2\sigma^2} \right),$$

or as a Julia function:
"""

# ╔═╡ 48530f0d-49b4-4aec-8109-d69f1ef7f0ee
md"""
Write a function `gauss` that takes `σ` as a keyword argument and implements this function.
"""

# ╔═╡ beb62fda-38a6-4528-a176-cfb726f4b5bd
gauss(x::Real; σ=1) = 1 / sqrt(2π*σ^2) * exp(-x^2 / (2 * σ^2))

# ╔═╡ f0d55cec-2e81-4cbb-b166-2cf4f2a0f43f
md"""
We need to **sample** (i.e. evaluate) this at each pixel in an interval of length $2n+1$,
and then **normalize** so that the sum of the resulting kernel is 1.
"""

# ╔═╡ f8bd22b8-ee14-11ea-04aa-ab16fd01826e
md"""
You can edit the cell above to test your kernel function!

Let's try applying it in a convolution.
"""

# ╔═╡ 2a9dd06a-ee13-11ea-3f84-67bb309c77a8
@bind gaussian_kernel_size_1D Slider(0:6)

# ╔═╡ 8f70cbfc-80f6-11eb-078c-337a02409c6b
x = [0.343, 0.23121, 1.0, 1.0]

# ╔═╡ 6f74c790-80f6-11eb-3200-717715b41e3d
colored_line(x)

# ╔═╡ ce24e486-df27-4780-bc57-d3bf7bee83bb
function create_bar()
	x = zeros(100)
	x[41:60] .= 1
	x
end

# ╔═╡ b01858b6-edf3-11ea-0826-938d33c19a43
md"""

## **Exercise 2** - _Convolutions in 2D_

Now let's move to 2D images. The convolution is then given by a **kernel matrix** $K$:

$$M'_{i, j} = \sum_{k, l}  \, M_{i- k, j - l} \, K_{k, l},$$

where the sum is over the possible values of $k$ and $l$ in the window. Again we think of the window as being *centered* at $(i, j)$.

A common notation for this operation is $\star$:

```math
M' = M \star K
```
"""

# ╔═╡ 7c1bc062-ee15-11ea-30b1-1b1e76520f13
md"""
#### Exercise 2.1
👉 Write a new method for `extend` that takes a matrix `M` and indices `i` and `j`, and returns the closest element of the matrix.
"""

# ╔═╡ 7c2ec6c6-ee15-11ea-2d7d-0d9401a5e5d1
function extend(M::AbstractMatrix, i, j)

	i = max(i, 1)
	i = min(i, size(M,1))

	j = max(j, 1)
	j = min(j, size(M,2))

	return M[i,j]
end

# ╔═╡ d87175b6-7db0-11eb-08a4-456a34fa137e
extend([1],0)

# ╔═╡ 803905b2-ee09-11ea-2d52-e77ff79693b0
extend([5,6,7], 1)

# ╔═╡ 80479d98-ee09-11ea-169e-d166eef65874
extend([5,6,7], -8)

# ╔═╡ 805691ce-ee09-11ea-053d-6d2e299ee123
extend([5,6,7], 10)

# ╔═╡ 45c4da9a-ee0f-11ea-2c5b-1f6704559137
if extend(v,1) === missing
	missing
else
	a = [extend(example_vector, i) for i in -1:length(example_vector)+2]
	colored_line(a)
end

# ╔═╡ 807e5662-ee09-11ea-3005-21fdcc36b023
function box_blur(v::AbstractArray, l)
	blurred_image = [mean(extend(v, j) for j in i-l:i+l) for i in 1:length(v)]
	return blurred_image
end

# ╔═╡ 4f08ebe8-b781-4a32-a218-5ecd8338561d
colored_line(box_blur(example_vector, 2))

# ╔═╡ 808deca8-ee09-11ea-0ee3-1586fa1ce282
let
	try
		test_v = rand(n)
		original = copy(test_v)
		box_blur(test_v, 5)
		if test_v != original
			md"""
			!!! danger "Oopsie!"
			    It looks like your function _modifies_ `v`. Can you write it without doing so? Maybe you can use `copy`.
			"""
		end
	catch
	end
end

# ╔═╡ 7d80a1ea-a0a9-41b2-9cfe-a334717ab2f4
colored_line(box_blur(v,l_box))

# ╔═╡ bbe1a562-8d97-4112-a88a-c45c260f574d
let
	result = box_blur(v, box_kernel_l)
	colored_line(result)
end

# ╔═╡ 28e20950-ee0c-11ea-0e0a-b5f2e570b56e
function convolve(v::AbstractVector, k)
	l = (length(k)-1)÷2
	offset_k = OffsetArray(k,-l:l)
	convolved_vector = [sum(offset_k[j]*extend(v, i+j) for j in -l:l) for i in 1:length(v)]
	return convolved_vector
end

# ╔═╡ 9afc4dca-ee16-11ea-354f-1d827aaa61d2
md"_Let's test it!_"

# ╔═╡ cf6b05e2-ee16-11ea-3317-8919565cb56e
small_image = Gray.(rand(5,5))

# ╔═╡ e3616062-ee27-11ea-04a9-b9ec60842a64
md"- Extended with `0`:"

# ╔═╡ e5b6cd34-ee27-11ea-0d60-bd4796540b18
[get(small_image, (i, j), Gray(0)) for (i,j) in Iterators.product(-1:7,-1:7)]

# ╔═╡ b4e98589-f221-4922-b11e-364d72d0788e

# ╔═╡ d06ea762-ee27-11ea-2e9c-1bcff86a3fe0
md"- Extended with your `extend` function:"

# ╔═╡ e1dc0622-ee16-11ea-274a-3b6ec9e15ab5
[extend(small_image, i, j) for (i,j) in Iterators.product(-1:7,-1:7)]

# ╔═╡ 4bbea325-35f8-4a51-bd66-153aba4aed96
md"""
### Extending Philip
"""

# ╔═╡ c4f5a867-74ba-4106-91d4-195f6ae644d0
url = "https://user-images.githubusercontent.com/6933510/107239146-dcc3fd00-6a28-11eb-8c7b-41aaf6618935.png"

# ╔═╡ c825ebe2-511b-43ba-afdf-6226dbac48d2
philip_filename = download(url) # download to a local file. The filename is returned

# ╔═╡ 2701ab0c-b91d-47fe-b36b-4e0036ecd4aa
philip = load(philip_filename);

# ╔═╡ 84a48984-9adb-40ab-a1f1-1ab7b76c9a19
philip_head = philip[470:800, 140:410];

# ╔═╡ 3cd535e4-ee26-11ea-2482-fb4ad43dda19
[
	extend(philip_head, i, j) for
		i in -50:size(philip_head,1)+51,
		j in -50:size(philip_head,2)+51
]

# ╔═╡ 7c41f0ca-ee15-11ea-05fb-d97a836659af
md"""
#### Exercise 2.2
👉 Implement a new method `convolve(M, K)` that applies a convolution to a 2D array `M`, using a 2D kernel `K`. Use your new method `extend` from the last exercise.
"""

# ╔═╡ 8b96e0bc-ee15-11ea-11cd-cfecea7075a0
function convolve(M::AbstractMatrix, K::AbstractMatrix)
	rows, cols = size(K)
	r = (rows-1)÷2
	c = (cols-1)÷2
	offset_kernel = OffsetArray(K, -r:r, -c:c)
	convolved_matrix = [sum([extend(M, i+m, j+n) * offset_kernel[i,j] for i in -r:r, j in -c:c]) for m in 1:size(M,1), n in 1:size(M,2)]

	return convolved_matrix
end

# ╔═╡ 93284f92-ee12-11ea-0342-833b1a30625c
test_convolution = let
	v = [1, 10, 100, 1000, 10000]
	k = [1, 1, 0]
	convolve(v, k)
end

# ╔═╡ 5eea882c-ee13-11ea-0d56-af81ecd30a4a
colored_line(test_convolution)

# ╔═╡ 338b1c3f-f071-4f80-86c0-a82c17349828
let
	result = convolve(v, box_blur_kernel_test)
	colored_line(result)
end

# ╔═╡ 5a5135c6-ee1e-11ea-05dc-eb0c683c2ce5
md"_Let's test it out! 🎃_"

# ╔═╡ 577c6daa-ee1e-11ea-1275-b7abc7a27d73
test_image_with_border = [get(small_image, (i, j), Gray(0)) for (i,j) in Iterators.product(-1:7,-1:7)]

# ╔═╡ 275a99c8-ee1e-11ea-0a76-93e3618c9588
K_test = [
	0   0  0
	1/2 0  1/2
	0   0  0
]

# ╔═╡ 42dfa206-ee1e-11ea-1fcd-21671042064c
new_image = convolve(test_image_with_border, K_test)

# ╔═╡ 2ac28130-810c-11eb-2cb1-117579142578
typeof(new_image)

# ╔═╡ 1d3b604a-810c-11eb-27fe-538cba860476
typeof(test_image_with_border)

# ╔═╡ 9ce383c8-810b-11eb-0c4c-fd184b9f6e3a
typeof(philip_head)

# ╔═╡ 4739eae4-810a-11eb-07f8-0b47904d9efe
convolve(philip_head, K_test)

# ╔═╡ 6e53c2e6-ee1e-11ea-21bd-c9c05381be07
md"_Edit_ `K_test` _to create your own test case!_"

# ╔═╡ e7f8b41a-ee25-11ea-287a-e75d33fbd98b
[philip_head convolve(philip_head, K_test)]

# ╔═╡ 8a335044-ee19-11ea-0255-b9391246d231
md"""
---

You can create all sorts of effects by choosing the kernel in a smart way. Today, we will implement two special kernels, to produce a **Gaussian blur** and a **Sobel edge detection** filter.

Make sure that you have watched the lecture about convolutions!
"""

# ╔═╡ 79eb0775-3582-446b-996a-0b64301394d0
md"""
#### Exercise 2.3
The 2D Gaussian kernel will be defined using

$$G(x,y)=\frac{1}{2\pi \sigma^2}\exp\left(\frac{-(x^2+y^2)}{2\sigma^2}\right)$$

How can you express this mathematically using the 1D Gaussian function that we defined before?
"""

# ╔═╡ f4d9fd6f-0f1b-4dec-ae68-e61550cee790
gauss(x, y; σ=1) = 2π*σ^2 * gauss(x; σ=σ) * gauss(y; σ=σ)

# ╔═╡ 1c8b4658-ee0c-11ea-2ede-9b9ed7d3125e
function gaussian_kernel_1D(n; σ = 1)
	gaussian_kernel = [gauss(x,σ) for x in -n:n]
	return gaussian_kernel./sum(gaussian_kernel)
end

# ╔═╡ a6149507-d5ba-45c1-896a-3487070d36ec
colored_line(gaussian_kernel_1D(4; σ=5))

# ╔═╡ 38eb92f6-ee13-11ea-14d7-a503ac04302e
test_gauss_1D_a = let
	k = gaussian_kernel_1D(gaussian_kernel_size_1D)

	if k !== missing
		convolve(x, k)
	end
end

# ╔═╡ b424e2aa-ee14-11ea-33fa-35491e0b9c9d
colored_line(test_gauss_1D_a)

# ╔═╡ 24c21c7c-ee14-11ea-1512-677980db1288
test_gauss_1D_b = let
	v = create_bar()
	k = gaussian_kernel_1D(gaussian_kernel_size_1D)

	if k !== missing
		convolve(v, k)
	end
end

# ╔═╡ bc1c20a4-ee14-11ea-3525-63c9fa78f089
colored_line(test_gauss_1D_b)

# ╔═╡ 7c50ea80-ee15-11ea-328f-6b4e4ff20b7e
md"""
👉 Write a function that applies a **Gaussian blur** to an image. Use your previous functions, and add cells to write helper functions as needed!
"""

# ╔═╡ 48ebad62-8107-11eb-34e9-09890a16e6ea
function gaussian_kernel_2D(n, σ)
	gaussian_kernel = [gauss(x, y, σ=σ) for x in -n:n, y in -n:n]
	return gaussian_kernel./sum(gaussian_kernel)
end

# ╔═╡ aad67fd0-ee15-11ea-00d4-274ec3cda3a3
function with_gaussian_blur(image; σ=3, l=5)
	return convolve(image, gaussian_kernel_2D(l, σ))
end

# ╔═╡ 8ae59674-ee18-11ea-3815-f50713d0fa08
md"_Let's make it interactive. 💫_"

# ╔═╡ 96146b16-79ea-401f-b8ba-e05663a18bd8
@bind face_σ Slider(0.1:0.1:10; show_value=true)

# ╔═╡ 2cc745ce-e145-4428-af3b-926fba271b67
@bind face_l Slider(0:20; show_value=true)

# ╔═╡ a75701c4-ee18-11ea-2863-d3042e71a68b
gaussian_blurred_image = with_gaussian_blur(philip_head; σ=face_σ, l=face_l)

# ╔═╡ 9cb8c8a0-810f-11eb-375b-7b4de6432ebd
[philip_head gaussian_blurred_image]

# ╔═╡ 364a0a42-8110-11eb-3907-e9cfc000a302
typeof(gaussian_blurred_image)

# ╔═╡ d5ffc6ab-156b-4d43-ac3d-1947d0176e7f
md"""
When you set `face_σ` to a low number (e.g. `2.0`), what effect does `face_l` have? And vice versa?
"""

# ╔═╡ 7c6642a6-ee15-11ea-0526-a1aac4286cdd
md"""
#### Exercise 2.4
👉 Create a **Sobel edge detection filter**.

Here, we will need to create two filters that separately detect edges in the horizontal and vertical directions, given by the following kernels:

```math
G_x = \begin{bmatrix}
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1 \\
\end{bmatrix};
\qquad
G_y = \begin{bmatrix}
1 & 2 & 1 \\
0 & 0 & 0 \\
-1 & -2 & -1 \\
\end{bmatrix}
```

We can think of these filterrs as derivatives in the $x$ and $y$ directions, as we discussed in lectures.

Then we combine them by finding the magnitude of the **gradient** (in the sense of multivariate calculus) by defining

$$G_\text{total} = \sqrt{G_x^2 + G_y^2},$$

where each operation applies *element-wise* on the matrices.

Use your previous functions, and add cells to write helper functions as needed!
"""

# ╔═╡ 9767e240-8110-11eb-3869-cbfed5e5918d
function horizontal_edge_detect_kernel()
	return [-0.125  -0.25  -0.125;
  0.0     0.0    0.0;
  0.125   0.25   0.125]
end

# ╔═╡ fe8eb9d0-8110-11eb-0fe6-d91ef91f892e
function vertical_edge_detect_kernel()
	return [1 2 1 ; 0 0 0 ; -1 -2 -1]
end

# ╔═╡ 9eeb876c-ee15-11ea-1794-d3ea79f47b75
function with_sobel_edge_detect(image)
	kernel_x = horizontal_edge_detect_kernel()
	kernel_y = vertical_edge_detect_kernel()

	x_edges = convolve(image, kernel_x)
	y_edges = convolve(image,kernel_y)

	y_edges

	edges_sqrt = sqrt.(x_edges.^2 + y_edges.^2)
	edges = x_edges + y_edges
	[edges_sqrt edges]
	# typeof(edges)
	# return edges
end

# ╔═╡ 525e8402-818c-11eb-1f8b-05cd978cffb3
Gray.(with_sobel_edge_detect(Gray.(philip_head)))

# ╔═╡ 260dc208-8113-11eb-05c9-33108c830d89
[Gray.(philip_head) Gray.(with_sobel_edge_detect(Gray.(philip_head)))]

# ╔═╡ d2063f66-818b-11eb-3762-e734ed57af8f
sobel_hor = Kernel.sobel()[1]

# ╔═╡ 10b2e046-818c-11eb-0f8f-af52231e7fca
convolve(philip_head, sobel_hor)

# ╔═╡ f864d0b4-818b-11eb-2c48-7b334976ae22
typeof(horizontal_edge_detect_kernel())

# ╔═╡ 8ffe16ce-ee20-11ea-18bd-15640f94b839
if student.kerberos_id === "jazz"
	md"""
!!! danger "Oops!"
    **Before you submit**, remember to fill in your name and kerberos ID at the top of this notebook!
	"""
end

# ╔═╡ 2d9f3ae4-9e4c-49ce-aab0-5f87aba85c3e
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 5516c800-edee-11ea-12cf-3f8c082ef0ef
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ ea435e58-ee11-11ea-3785-01af8dd72360
hint(md"Have a look at the lecture notes to see examples of adding interactivity with a slider. You can read the Interactivity and the PlutoUI sample notebooks to learn more, you can find them in Pluto's main menu. _(Right click the Pluto logo in the top left -> Open in new tab)_.")

# ╔═╡ 32a07f1d-93cd-4bf3-bac1-91afa6bb88a6
md"""
You can use the `÷` operator (you type `\div<TAB>` to get it with autocomplete) to do _integer division_. For example:

```julia
8 / 6 ≈ 1.3333333 # a floating point number!

8 // 6 == 4 // 3  # a fraction!

8 ÷ 6 == 1        # an integer!
```
""" |> hint

# ╔═╡ 649df270-ee24-11ea-397e-79c4355e38db
hint(md"`num_rows, num_columns = size(M)`")

# ╔═╡ 0cabed84-ee1e-11ea-11c1-7d8a4b4ad1af
hint(md"`num_rows, num_columns = size(K)`")

# ╔═╡ 9def5f32-ee15-11ea-1f74-f7e6690f2efa
hint(md"Can we just copy the 1D code? What is different in 2D?")

# ╔═╡ 57360a7a-edee-11ea-0c28-91463ece500d
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ dcb8324c-edee-11ea-17ff-375ff5078f43
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ 58af703c-edee-11ea-2963-f52e78fc2412
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ f3d00a9a-edf3-11ea-07b3-1db5c6d0b3cf
yays = [md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next exercise."]

# ╔═╡ 5aa9dfb2-edee-11ea-3754-c368fb40637c
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ f0c3e99d-9eb9-459e-917a-c2338af6683c
let
	result = gaussian_kernel_1D(5)

	if ismissing(result)
		still_missing()
	elseif isnothing(result)
		keep_working(md"Did you forget to write `return`?")
	elseif !(result isa AbstractVector)
		keep_working(md"The returned object is not a `Vector`.")
	elseif size(result) != (11,)
		hint(md"The returned vector has the wrong dimensions.")
	elseif !(sum(result) ≈ 1.0)
		keep_working(md"You need to _normalize_ the result.")
	elseif gaussian_kernel_1D(3; σ=1) == gaussian_kernel_1D(3; σ=2)
		keep_working(md"Use the keyword argument `σ` in your function.")
	else
		correct()
	end
end

# ╔═╡ 74d44e22-edee-11ea-09a0-69aa0aba3281
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ bcf98dfc-ee1b-11ea-21d0-c14439500971
if !@isdefined(extend)
	not_defined(:extend)
else
	let
		result = extend([6,7],-10)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != 6 || extend([6,7],10) != 7
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 7ffd14f8-ee1d-11ea-0343-b54fb0333aea
if !@isdefined(convolve)
	not_defined(:convolve)
else
	let
		x = [1, 10, 100]
		result = convolve(x, [0, 1, 1])
		shouldbe = [11, 110, 200]
		shouldbe2 = [2, 11, 110]

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa AbstractVector)
			keep_working(md"The returned object is not a `Vector`.")
		elseif size(result) != size(x)
			keep_working(md"The returned vector has the wrong dimensions.")
		elseif result != shouldbe && result != shouldbe2
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ d93fa3f6-c361-4dfd-a2ea-f38e682bcd6a
if !@isdefined(box_blur_kernel)
	not_defined(:box_blur_kernel)
else
	let
		result = box_blur_kernel(2)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa AbstractVector)
			keep_working(md"The returned object is not a `Vector`.")
		elseif size(result) != (5,)
			hint(md"The returned vector has the wrong dimensions.")
		else

			x = [1, 10, 100]
			result1 = box_blur(x, 2)
			result2 = convolve(x, result)

			if result1 ≈ result2
				correct()
			else
				keep_working()
			end
		end
	end
end

# ╔═╡ efd1ceb4-ee1c-11ea-350e-f7e3ea059024
if !@isdefined(extend)
	not_defined(:extend)
else
	let
		input = [42 37; 1 0]
		result = extend(input, -2, -2)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != 42 || extend(input, -1, 3) != 37
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 115ded8c-ee0a-11ea-3493-89487315feb7
bigbreak = html"<br><br><br><br><br>";

# ╔═╡ 54056a02-ee0a-11ea-101f-47feb6623bec
bigbreak

# ╔═╡ a3067222-a83a-47b8-91c3-24ad78dd65c5
bigbreak

# ╔═╡ 27847dc4-ee0a-11ea-0651-ebbbb3cfd58c
bigbreak

# ╔═╡ 0001f782-ee0e-11ea-1fb4-2b5ef3d241e2
bigbreak

# ╔═╡ 5842895a-ee10-11ea-119d-81e4c4c8c53b
bigbreak

# ╔═╡ dfb7c6be-ee0d-11ea-194e-9758857f7b20
function camera_input(;max_size=200, default_url="https://i.imgur.com/SUmi94P.png")
"""
<span class="pl-image waiting-for-permission">
<style>

	.pl-image.popped-out {
		position: fixed;
		top: 0;
		right: 0;
		z-index: 5;
	}

	.pl-image #video-container {
		width: 250px;
	}

	.pl-image video {
		border-radius: 1rem 1rem 0 0;
	}
	.pl-image.waiting-for-permission #video-container {
		display: none;
	}
	.pl-image #prompt {
		display: none;
	}
	.pl-image.waiting-for-permission #prompt {
		width: 250px;
		height: 200px;
		display: grid;
		place-items: center;
		font-family: monospace;
		font-weight: bold;
		text-decoration: underline;
		cursor: pointer;
		border: 5px dashed rgba(0,0,0,.5);
	}

	.pl-image video {
		display: block;
	}
	.pl-image .bar {
		width: inherit;
		display: flex;
		z-index: 6;
	}
	.pl-image .bar#top {
		position: absolute;
		flex-direction: column;
	}

	.pl-image .bar#bottom {
		background: black;
		border-radius: 0 0 1rem 1rem;
	}
	.pl-image .bar button {
		flex: 0 0 auto;
		background: rgba(255,255,255,.8);
		border: none;
		width: 2rem;
		height: 2rem;
		border-radius: 100%;
		cursor: pointer;
		z-index: 7;
	}
	.pl-image .bar button#shutter {
		width: 3rem;
		height: 3rem;
		margin: -1.5rem auto .2rem auto;
	}

	.pl-image video.takepicture {
		animation: pictureflash 200ms linear;
	}

	@keyframes pictureflash {
		0% {
			filter: grayscale(1.0) contrast(2.0);
		}

		100% {
			filter: grayscale(0.0) contrast(1.0);
		}
	}
</style>

	<div id="video-container">
		<div id="top" class="bar">
			<button id="stop" title="Stop video">✖</button>
			<button id="pop-out" title="Pop out/pop in">⏏</button>
		</div>
		<video playsinline autoplay></video>
		<div id="bottom" class="bar">
		<button id="shutter" title="Click to take a picture">📷</button>
		</div>
	</div>

	<div id="prompt">
		<span>
		Enable webcam
		</span>
	</div>

<script>
	// based on https://github.com/fonsp/printi-static (by the same author)

	const span = currentScript.parentElement
	const video = span.querySelector("video")
	const popout = span.querySelector("button#pop-out")
	const stop = span.querySelector("button#stop")
	const shutter = span.querySelector("button#shutter")
	const prompt = span.querySelector(".pl-image #prompt")

	const maxsize = $(max_size)

	const send_source = (source, src_width, src_height) => {
		const scale = Math.min(1.0, maxsize / src_width, maxsize / src_height)

		const width = Math.floor(src_width * scale)
		const height = Math.floor(src_height * scale)

		const canvas = html`<canvas width=\${width} height=\${height}>`
		const ctx = canvas.getContext("2d")
		ctx.drawImage(source, 0, 0, width, height)

		span.value = {
			width: width,
			height: height,
			data: ctx.getImageData(0, 0, width, height).data,
		}
		span.dispatchEvent(new CustomEvent("input"))
	}

	const clear_camera = () => {
		window.stream.getTracks().forEach(s => s.stop());
		video.srcObject = null;

		span.classList.add("waiting-for-permission");
	}

	prompt.onclick = () => {
		navigator.mediaDevices.getUserMedia({
			audio: false,
			video: {
				facingMode: "environment",
			},
		}).then(function(stream) {

			stream.onend = console.log

			window.stream = stream
			video.srcObject = stream
			window.cameraConnected = true
			video.controls = false
			video.play()
			video.controls = false

			span.classList.remove("waiting-for-permission");

		}).catch(function(error) {
			console.log(error)
		});
	}
	stop.onclick = () => {
		clear_camera()
	}
	popout.onclick = () => {
		span.classList.toggle("popped-out")
	}

	shutter.onclick = () => {
		const cl = video.classList
		cl.remove("takepicture")
		void video.offsetHeight
		cl.add("takepicture")
		video.play()
		video.controls = false
		console.log(video)
		send_source(video, video.videoWidth, video.videoHeight)
	}

	document.addEventListener("visibilitychange", () => {
		if (document.visibilityState != "visible") {
			clear_camera()
		}
	})

	// Set a default image

	const img = html`<img crossOrigin="anonymous">`

	img.onload = () => {
	console.log("helloo")
		send_source(img, img.width, img.height)
	}
	img.src = "$(default_url)"
	console.log(img)
</script>
</span>
""" |> HTML
end

# ╔═╡ 94c0798e-ee18-11ea-3212-1533753eabb6
@bind gauss_raw_camera_data camera_input(;max_size=100)

# ╔═╡ e7b91064-8112-11eb-06c1-cd57fd57133f
gauss_raw_camera_data

# ╔═╡ 1a0324de-ee19-11ea-1d4d-db37f4136ad3
@bind sobel_raw_camera_data camera_input(;max_size=200)

# ╔═╡ e15ad330-ee0d-11ea-25b6-1b1b3f3d7888

function process_raw_camera_data(raw_camera_data)
	# the raw image data is a long byte array, we need to transform it into something
	# more "Julian" - something with more _structure_.

	# The encoding of the raw byte stream is:
	# every 4 bytes is a single pixel
	# every pixel has 4 values: Red, Green, Blue, Alpha
	# (we ignore alpha for this notebook)

	# So to get the red values for each pixel, we take every 4th value, starting at
	# the 1st:
	reds_flat = UInt8.(raw_camera_data["data"][1:4:end])
	greens_flat = UInt8.(raw_camera_data["data"][2:4:end])
	blues_flat = UInt8.(raw_camera_data["data"][3:4:end])

	# but these are still 1-dimensional arrays, nicknamed 'flat' arrays
	# We will 'reshape' this into 2D arrays:

	width = raw_camera_data["width"]
	height = raw_camera_data["height"]

	# shuffle and flip to get it in the right shape
	reds = reshape(reds_flat, (width, height))' / 255.0
	greens = reshape(greens_flat, (width, height))' / 255.0
	blues = reshape(blues_flat, (width, height))' / 255.0

	# we have our 2D array for each color
	# Let's create a single 2D array, where each value contains the R, G and B value of
	# that pixel

	RGB.(reds, greens, blues)
end

# ╔═╡ f461f5f2-ee18-11ea-3d03-95f57f9bf09e
gauss_camera_image = process_raw_camera_data(gauss_raw_camera_data);

# ╔═╡ 1ff6b5cc-ee19-11ea-2ca8-7f00c204f587
sobel_camera_image = Gray.(process_raw_camera_data(sobel_raw_camera_data));

# ╔═╡ 1bf94c00-ee19-11ea-0e3c-e12bc68d8e28
[sobel_camera_image Gray.(with_sobel_edge_detect(Gray.(sobel_camera_image)))]

# ╔═╡ Cell order:
# ╟─83eb9ca0-ed68-11ea-0bc5-99a09c68f867
# ╟─8ef13896-ed68-11ea-160b-3550eeabbd7d
# ╟─ac8ff080-ed61-11ea-3650-d9df06123e1f
# ╠═911ccbce-ed68-11ea-3606-0384e7580d7c
# ╟─5f95e01a-ee0a-11ea-030c-9dba276aba92
# ╠═65780f00-ed6b-11ea-1ecf-8b35523a7ac0
# ╟─f7a6d7c3-37b9-437d-8b8e-853665ddbae3
# ╟─54056a02-ee0a-11ea-101f-47feb6623bec
# ╟─e08781fa-ed61-11ea-13ae-91a49b5eb74a
# ╟─a3067222-a83a-47b8-91c3-24ad78dd65c5
# ╟─80108d80-ee09-11ea-0368-31546eb0d3cc
# ╠═7fcd6230-ee09-11ea-314f-a542d00d582e
# ╠═7fdb34dc-ee09-11ea-366b-ffe10d1aa845
# ╟─7fe9153e-ee09-11ea-15b3-6f24fcc20734
# ╠═01070e28-ee0f-11ea-1928-a7919d452bdd
# ╠═ff70782e-e8d2-4281-9b24-d45c925f55e2
# ╟─7522f81e-ee1c-11ea-35af-a17eb257ff1a
# ╟─801d90c0-ee09-11ea-28d6-61b806de26dc
# ╠═802bec56-ee09-11ea-043e-51cf1db02a34
# ╠═d87175b6-7db0-11eb-08a4-456a34fa137e
# ╟─b7f3994c-ee1b-11ea-211a-d144db8eafc2
# ╠═803905b2-ee09-11ea-2d52-e77ff79693b0
# ╠═80479d98-ee09-11ea-169e-d166eef65874
# ╠═805691ce-ee09-11ea-053d-6d2e299ee123
# ╟─bcf98dfc-ee1b-11ea-21d0-c14439500971
# ╠═3492b164-7065-48e8-978b-6c96b965d376
# ╠═02123165-2a0a-49a8-b7a9-458955523511
# ╟─806e5766-ee0f-11ea-1efc-d753cd83d086
# ╠═38da843a-ee0f-11ea-01df-bfa8b1317d36
# ╟─9bde9f92-ee0f-11ea-27f8-ffef5fce2b3c
# ╠═45c4da9a-ee0f-11ea-2c5b-1f6704559137
# ╟─431ba330-0f72-416a-92e9-55f51ff3bcd1
# ╠═5fdc5d0d-a52c-476e-b3b5-3b6364b706e4
# ╟─e84c9cc2-e6e1-46f1-bf4e-9605da5e6f4a
# ╠═807e5662-ee09-11ea-3005-21fdcc36b023
# ╠═ca91bf1a-80ed-11eb-3be0-cd60c9adae9d
# ╠═4f08ebe8-b781-4a32-a218-5ecd8338561d
# ╟─808deca8-ee09-11ea-0ee3-1586fa1ce282
# ╟─809f5330-ee09-11ea-0e5b-415044b6ac1f
# ╟─e555a7e6-f11a-43ac-8218-6d832f0ce251
# ╠═302f0842-453f-47bd-a74c-7942d8c96485
# ╠═7d80a1ea-a0a9-41b2-9cfe-a334717ab2f4
# ╟─ea435e58-ee11-11ea-3785-01af8dd72360
# ╟─80ab64f4-ee09-11ea-29b4-498112ed0799
# ╠═28e20950-ee0c-11ea-0e0a-b5f2e570b56e
# ╟─32a07f1d-93cd-4bf3-bac1-91afa6bb88a6
# ╟─5eea882c-ee13-11ea-0d56-af81ecd30a4a
# ╠═93284f92-ee12-11ea-0342-833b1a30625c
# ╟─cf73f9f8-ee12-11ea-39ae-0107e9107ef5
# ╟─7ffd14f8-ee1d-11ea-0343-b54fb0333aea
# ╟─fa463b71-5aa4-44a3-a67b-6b0776236243
# ╠═8a7d3cfd-6f19-43f0-ae16-d5a236f148e7
# ╟─a34d1ad8-3776-4bc4-93e5-72cfffc54f15
# ╠═971a801d-9c46-417a-ad31-1144894fb4e1
# ╟─5f13b1a5-8c7d-47c9-b96a-a09faf38fe5e
# ╠═338b1c3f-f071-4f80-86c0-a82c17349828
# ╠═bbe1a562-8d97-4112-a88a-c45c260f574d
# ╟─d93fa3f6-c361-4dfd-a2ea-f38e682bcd6a
# ╟─03f91a22-1c3e-4c42-9d78-1ee36851a120
# ╟─48530f0d-49b4-4aec-8109-d69f1ef7f0ee
# ╠═beb62fda-38a6-4528-a176-cfb726f4b5bd
# ╟─f0d55cec-2e81-4cbb-b166-2cf4f2a0f43f
# ╠═1c8b4658-ee0c-11ea-2ede-9b9ed7d3125e
# ╟─f0c3e99d-9eb9-459e-917a-c2338af6683c
# ╠═a6149507-d5ba-45c1-896a-3487070d36ec
# ╟─f8bd22b8-ee14-11ea-04aa-ab16fd01826e
# ╠═2a9dd06a-ee13-11ea-3f84-67bb309c77a8
# ╠═8f70cbfc-80f6-11eb-078c-337a02409c6b
# ╠═6f74c790-80f6-11eb-3200-717715b41e3d
# ╠═b424e2aa-ee14-11ea-33fa-35491e0b9c9d
# ╠═38eb92f6-ee13-11ea-14d7-a503ac04302e
# ╠═bc1c20a4-ee14-11ea-3525-63c9fa78f089
# ╠═24c21c7c-ee14-11ea-1512-677980db1288
# ╟─ce24e486-df27-4780-bc57-d3bf7bee83bb
# ╟─27847dc4-ee0a-11ea-0651-ebbbb3cfd58c
# ╟─b01858b6-edf3-11ea-0826-938d33c19a43
# ╟─7c1bc062-ee15-11ea-30b1-1b1e76520f13
# ╠═7c2ec6c6-ee15-11ea-2d7d-0d9401a5e5d1
# ╟─649df270-ee24-11ea-397e-79c4355e38db
# ╟─9afc4dca-ee16-11ea-354f-1d827aaa61d2
# ╠═cf6b05e2-ee16-11ea-3317-8919565cb56e
# ╟─e3616062-ee27-11ea-04a9-b9ec60842a64
# ╟─e5b6cd34-ee27-11ea-0d60-bd4796540b18
# ╟─b4e98589-f221-4922-b11e-364d72d0788e
# ╟─d06ea762-ee27-11ea-2e9c-1bcff86a3fe0
# ╟─e1dc0622-ee16-11ea-274a-3b6ec9e15ab5
# ╟─efd1ceb4-ee1c-11ea-350e-f7e3ea059024
# ╟─4bbea325-35f8-4a51-bd66-153aba4aed96
# ╠═c4f5a867-74ba-4106-91d4-195f6ae644d0
# ╠═c825ebe2-511b-43ba-afdf-6226dbac48d2
# ╠═2701ab0c-b91d-47fe-b36b-4e0036ecd4aa
# ╠═84a48984-9adb-40ab-a1f1-1ab7b76c9a19
# ╠═3cd535e4-ee26-11ea-2482-fb4ad43dda19
# ╟─7c41f0ca-ee15-11ea-05fb-d97a836659af
# ╠═8b96e0bc-ee15-11ea-11cd-cfecea7075a0
# ╟─0cabed84-ee1e-11ea-11c1-7d8a4b4ad1af
# ╟─5a5135c6-ee1e-11ea-05dc-eb0c683c2ce5
# ╟─577c6daa-ee1e-11ea-1275-b7abc7a27d73
# ╠═275a99c8-ee1e-11ea-0a76-93e3618c9588
# ╠═42dfa206-ee1e-11ea-1fcd-21671042064c
# ╠═2ac28130-810c-11eb-2cb1-117579142578
# ╠═1d3b604a-810c-11eb-27fe-538cba860476
# ╠═9ce383c8-810b-11eb-0c4c-fd184b9f6e3a
# ╠═4739eae4-810a-11eb-07f8-0b47904d9efe
# ╟─6e53c2e6-ee1e-11ea-21bd-c9c05381be07
# ╠═e7f8b41a-ee25-11ea-287a-e75d33fbd98b
# ╟─8a335044-ee19-11ea-0255-b9391246d231
# ╟─79eb0775-3582-446b-996a-0b64301394d0
# ╠═f4d9fd6f-0f1b-4dec-ae68-e61550cee790
# ╟─7c50ea80-ee15-11ea-328f-6b4e4ff20b7e
# ╠═48ebad62-8107-11eb-34e9-09890a16e6ea
# ╠═aad67fd0-ee15-11ea-00d4-274ec3cda3a3
# ╟─9def5f32-ee15-11ea-1f74-f7e6690f2efa
# ╟─8ae59674-ee18-11ea-3815-f50713d0fa08
# ╠═94c0798e-ee18-11ea-3212-1533753eabb6
# ╠═e7b91064-8112-11eb-06c1-cd57fd57133f
# ╠═a75701c4-ee18-11ea-2863-d3042e71a68b
# ╠═9cb8c8a0-810f-11eb-375b-7b4de6432ebd
# ╠═364a0a42-8110-11eb-3907-e9cfc000a302
# ╠═96146b16-79ea-401f-b8ba-e05663a18bd8
# ╠═2cc745ce-e145-4428-af3b-926fba271b67
# ╟─d5ffc6ab-156b-4d43-ac3d-1947d0176e7f
# ╟─f461f5f2-ee18-11ea-3d03-95f57f9bf09e
# ╟─7c6642a6-ee15-11ea-0526-a1aac4286cdd
# ╠═9767e240-8110-11eb-3869-cbfed5e5918d
# ╠═fe8eb9d0-8110-11eb-0fe6-d91ef91f892e
# ╠═9eeb876c-ee15-11ea-1794-d3ea79f47b75
# ╠═525e8402-818c-11eb-1f8b-05cd978cffb3
# ╠═1a0324de-ee19-11ea-1d4d-db37f4136ad3
# ╠═1bf94c00-ee19-11ea-0e3c-e12bc68d8e28
# ╠═260dc208-8113-11eb-05c9-33108c830d89
# ╠═d2063f66-818b-11eb-3762-e734ed57af8f
# ╠═10b2e046-818c-11eb-0f8f-af52231e7fca
# ╠═f864d0b4-818b-11eb-2c48-7b334976ae22
# ╟─1ff6b5cc-ee19-11ea-2ca8-7f00c204f587
# ╟─0001f782-ee0e-11ea-1fb4-2b5ef3d241e2
# ╟─8ffe16ce-ee20-11ea-18bd-15640f94b839
# ╟─5842895a-ee10-11ea-119d-81e4c4c8c53b
# ╟─2d9f3ae4-9e4c-49ce-aab0-5f87aba85c3e
# ╟─5516c800-edee-11ea-12cf-3f8c082ef0ef
# ╟─57360a7a-edee-11ea-0c28-91463ece500d
# ╟─dcb8324c-edee-11ea-17ff-375ff5078f43
# ╟─58af703c-edee-11ea-2963-f52e78fc2412
# ╟─f3d00a9a-edf3-11ea-07b3-1db5c6d0b3cf
# ╟─5aa9dfb2-edee-11ea-3754-c368fb40637c
# ╟─74d44e22-edee-11ea-09a0-69aa0aba3281
# ╟─115ded8c-ee0a-11ea-3493-89487315feb7
# ╟─dfb7c6be-ee0d-11ea-194e-9758857f7b20
# ╟─e15ad330-ee0d-11ea-25b6-1b1b3f3d7888


text:
export ParamBitCircuit

import LogicCircuits: num_nodes, num_elements, num_features, num_leafs,
    nodes, elements, num_examples, isbinarydata, iscomplete

@inline num_features(x::BitMatrix) = size(x)[2]
@inline num_examples(x::BitMatrix) = size(x)[1]
@inline isbinarydata(x::BitMatrix) = true
@inline iscomplete(x::BitMatrix) = true

@inline num_classes(pbc::ParamBitCircuit, ::Type{LogisticCircuit}) = size(pbc.params)[2]

function ParamBitCircuit(lc::LogisticCircuit, data)
    nc = num_classes(lc)
    thetas::Vector{Vector{Float32}} = Vector{Vector{Float32}}()
    on_decision(n, cs, layer_id, decision_id, first_element, last_element) = begin
        if isnothing(n)
            # @assert first_element == last_element
            push!(thetas, zeros(Float32, nc))
            # println("here, some node is not part of the logistic circuit")
        else
            # @assert last_element - first_element + 1 == size(n.thetas, 1)
            # @assert size(n.thetas, 2) == nc
            for theta in eachrow(n.thetas)
                push!(thetas, theta)
            end
        end
    end
    bc = BitCircuit(lc, data; on_decision)
    ParamBitCircuit(bc, permutedims(hcat(thetas...), (2, 1)))
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
using BenchmarkTools
using Printf

# Get the dimension from the command line.
n, = size(ARGS)
if n < 1
    println("Usage: copyMatrix.jl dim")
    println("       --->      Please specify the dimension.")
    exit()
end

dimension = parse(Int, ARGS[1])

"""
    maxtrixcopy(A)

Perform the copy operations on matrix A with the given dimensions.
"""
function matrixcopy(A)
    N = size(A, 1)
    for j = 1:N, i = 1:N
        A[i, j, 1] = A[i, j, 2]
        A[i, j, 3] = A[i, j, 1]
        A[i, j, 2] = A[i, j, 3]
    end
end

println("-------------------------------")
println(@sprintf "Copy of matrix (loop) %d" dimension)
println("-------------------------------")

A = randn(dimension, dimension, 3)

@btime matrixcopy(A)

println(" ")

println("--------------------------")
println(@sprintf "Vectorized Copy of matrix %d" dimension)
println("--------------------------")

A = randn(dimension, dimension, 3)

@btime begin
    A[:, :, 1] = A[:, :, 2]
    A[:, :, 3] = A[:, :, 1]
    A[:, :, 2] = A[:, :, 3]
end

println(" ")


text:
"""
normalized water saturation
sws(sw::Real, sor::Real, swc::Real)
"""
function sws(sw::Real, sor::Real, swc::Real)
  if sw>swc && sw<1-sor
    res=(sw-swc)/(1-sor-swc)
  elseif sw>=1-sor
    res=1.0
  elseif sw<=swc
    res=0.0
  end
  return res
end

"""
Oil relative permeability functions for oil-water system
kro(sw::Real, kro0::Real, sor::Real, swc::Real, no::Real)
"""
function kro(sw::Real, kro0::Real, sor::Real, swc::Real, no::Real)
  if swc<=sw<=1-sor
    res=kro0*((1-sw-sor)/(1-sor-swc))^no
  elseif 0.0<sw<swc
    res=1+(kro0-1)/swc*sw
  elseif sw>1-sor
    res=0.0
  elseif sw<=0.0
    res=1.0
  end
  return res
end

"""
derivative of the oil relperm with respect to water saturation
dkrodsw(sw::Real, kro0::Real, sor::Real, swc::Real, no::Real)
"""
function dkrodsw(sw::Real, kro0::Real, sor::Real, swc::Real, no::Real)
  if swc<=sw<=1-sor
    res=-no*kro0/(1-sor-swc)*((1-sw-sor)/(1-sor-swc))^(no-1)
  elseif 0.0<sw<swc
    res=(kro0-1)/swc
  elseif sw>1-sor || sw<=0.0
    res=0.0
  end
  return res
end

"""
water rel-perm for oil-water system
krw(sw::Real, krw0::Real, sor::Real, swc::Real, nw::Real)
"""
function krw(sw::Real, krw0::Real, sor::Real, swc::Real, nw::Real)
  if swc<=sw<=1-sor
    res=krw0*((sw-swc)/(1-sor-swc))^nw
  elseif 1-sor<sw<1.0
    res=-(1-krw0)/sor*(1.0-sw)+1.0
  elseif sw<swc
    res=0.0
  elseif sw>=1.0
    res=1.0
  end
  return res
end

"""
derivative of the water relperm with respect to water saturation
dkrwdsw(sw::Real, krw0::Real, sor::Real, swc::Real, nw::Real)
"""
function dkrwdsw(sw::Real, krw0::Real, sor::Real, swc::Real, nw::Real)
  if swc<=sw<=1-sor
    res=nw*krw0/(1-sor-swc)*((sw-swc)/(1-sor-swc))^(nw-1)
  elseif 1-sor<sw<1.0
    res=(1-krw0)/sor
  elseif sw<swc || sw>=1.0
    res=0.0
  end
  return res
end

"""
capillary pressure curve for drainage
pc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e7)
"""
function pc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e7)
  sw0=swc+(1-labda*log(pc0/pce)+sqrt((-1+labda*log(pc0/pce))^2+4*swc/(1-swc)))/2*(1-swc)
  if sw>sw0
    res=pce*((sw-swc)/(1-swc))^(-1.0/labda)
  elseif 0.0<=sw<=sw0
    pcs=pce*((sw0-swc)/(1-swc))^(-1.0/labda)
    res=exp((log(pcs)-log(pc0))/sw0*(sw-sw0)+log(pcs))
  elseif sw<0
    res=pc0
  else
    res=0.0
  end
end

"""
derivative of the drainage capillary pressure with respect to water saturation
dpc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e7)
"""
function dpc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e7)
  sw0=swc+(1-labda*log(pc0/pce)+sqrt((-1+labda*log(pc0/pce))^2+4*swc/(1-swc)))/2*(1-swc)
  if sw>sw0
    res=-1.0/((1-swc)*labda)*pce*((sw-swc)/(1-swc))^(-1.0/labda-1)
  elseif 0.0<=sw<=sw0
    res=-1.0/((1-swc)*labda)*pce*((sw0-swc)/(1-swc))^(-1.0/labda-1)
  else
    res=0.0
  end
end

"""
second derivative of the drainage capillary pressure with respect to water saturation
dpc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e7)
"""
function d2pc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e7)
  sw0=swc+(1-labda*log(pc0/pce)+sqrt((-1+labda*log(pc0/pce))^2+4*swc/(1-swc)))/2*(1-swc)
  if sw>sw0
    res=-1.0/((1-swc)^2*labda)*(-1.0/labda-1)*pce*((sw-swc)/(1-swc))^(-1.0/labda-2.0)
  elseif 0.0<=sw<=sw0
    res=0.0
  else
    res=0.0
  end
end

"""
Imbibition capillary pressure curve
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e7, pc02::Real=1.0e7)
"""
function pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e7, pc02::Real=1.0e7)
  pc1=pc_drain(sw, pce, swc, labda=labda, pc0=pc01)
  pc2=pc_drain(1-sw, pce, sor, labda=labda, pc0=pc02)
  return (0.5*(1+cos(teta)))^b*pc1-(0.5*(1-cos(teta)))^b*pc2
end

"""
derivative of the imbibition capillary pressure curve w.r.t water saturation
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e7, pc02::Real=1.0e7)
"""
function dpc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda1::Real=2.4, labda2::Real=2.4, b::Real=0.6, pc01::Real=1.0e7, pc02::Real=1.0e7)
  dpc1=dpc_drain(sw, pce, swc, labda=labda1, pc0=pc01)
  dpc2=dpc_drain(1-sw, pce, sor, labda=labda2, pc0=pc02)
  return (0.5*(1+cos(teta)))^b*dpc1+(0.5*(1-cos(teta)))^b*dpc2
end

"""
second derivative of the imbibition capillary pressure curve w.r.t water saturation
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e7, pc02::Real=1.0e7)
"""
function d2pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda1::Real=2.4, labda2::Real=2.4, b::Real=0.6, pc01::Real=1.0e7, pc02::Real=1.0e7)
  dpc1=d2pc_drain(sw, pce, swc, labda=labda1, pc0=pc01)
  dpc2=d2pc_drain(1-sw, pce, sor, labda=labda2, pc0=pc02)
  return (0.5*(1+cos(teta)))^b*dpc1+(0.5*(1-cos(teta)))^b*dpc2
end

"""
SECOND approach (see the doc folder)
capillary pressure curve for drainage
pc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e6)
"""
function pc_drain2(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc_star::Real=1.0e6)
  sw_star=(1-swc)*exp(-labda*log(pc_star/pce))+swc
  pc0 = pc_star*exp(sw_star/(labda*(1-swc)))
  if sw>sw_star
    res=pce*((sw-swc)/(1.0-swc))^(-1.0/labda)
  elseif 0.0<=sw<=sw_star
    res=pc_star*exp((sw_star-sw)/(labda*(1.0-swc)))
  elseif sw<0
    res=pc0
  else
    res=0.0
  end
end

"""

derivative of the drainage capillary pressure with respect to water saturation
dpc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e6)
"""
function dpc_drain2(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc_star::Real=1.0e6)
  sw_star=(1-swc)*exp(-labda*log(pc_star/pce))+swc
  if sw>sw_star
    res=-1.0/((1-swc)*labda)*pce*((sw-swc)/(1-swc))^(-1.0/labda-1.0)
  elseif 0.0<=sw<=sw_star
    res=pc_star/(sw-sw_star)
  else
    res=0.0
  end
end

"""
Imbibition capillary pressure curve
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e6, pc02::Real=1.0e6)
"""
function pc_imb2(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc_star1::Real=1.0e6, pc_star2::Real=1.0e6)
  pc1=pc_drain2(sw, pce, swc, labda=labda, pc_star=pc_star1)
  pc2=pc_drain2(1-sw, pce, sor, labda=labda, pc_star=pc_star2)
  return (0.5*(1+cos(teta)))^b*pc1-(0.5*(1-cos(teta)))^b*pc2
end

"""
derivative of the imbibition capillary pressure curve w.r.t water saturation
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e6, pc02::Real=1.0e6)
"""
function dpc_imb2(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda1::Real=2.4, labda2::Real=2.4, b::Real=0.6, pc_star1::Real=1.0e6, pc_star2::Real=1.0e7)
  dpc1=dpc_drain2(sw, pce, swc, labda=labda1, pc_star=pc_star1)
  dpc2=dpc_drain2(1-sw, pce, sor, labda=labda2, pc_star=pc_star2)
  return (0.5*(1+cos(teta)))^b*dpc1+(0.5*(1-cos(teta)))^b*dpc2
end

"""
Third approach (see the doc folder)
capillary pressure curve for drainage
pc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e6)
"""
function pc_drain3(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc_star::Real=1.0e6, c::Real=1e-3)
  # pc_max  = pc_max_multiplyer*pc_star
  sw_star=(1-swc)*exp(-labda*log(pc_star/pce))+swc
  # m = -pce/(labda*(1-swc))*((sw_star-swc)/(1-swc))^(-1.0/labda-1.0)
  # c = 0.5*(-1.0+sqrt(1.0-(pc_max-pc_star)/(m*(exp(sw_star)-1))))
  b = -pce*c/(labda*(1-swc))*((sw_star-swc)/(1-swc))^(-1.0/labda-1.0)
  a = pc_star + b*exp(c)
  pc0 = a-b*exp(sw_star+c)
  if sw>sw_star
    res=pce*((sw-swc)/(1-swc))^(-1.0/labda)
  elseif 0.0<=sw<=sw_star
    res=a-b*exp(sw_star-sw+c)
  elseif sw<0
    res=pc0
  else
    res=0.0
  end
end

"""
derivative of the drainage capillary pressure with respect to water saturation
dpc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e6)
"""
function dpc_drain3(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc_star::Real=1.0e6, c::Real=1e-3)
  sw_star=(1-swc)*exp(-labda*log(pc_star/pce))+swc
  b = -pce*c/(labda*(1-swc))*((sw_star-swc)/(1-swc))^(-1.0/labda-1.0)
  a = pc_star + b*exp(c)
  pc0 = a-b*exp(sw_star+c)
  if sw>sw_star
    res=-1.0/((1-swc)*labda)*pce*((sw-swc)/(1-swc))^(-1.0/labda-1)
  elseif 0.0<=sw<=sw_star
    res=b/(sw_star-sw+c)
  else
    res=b/(sw_star+c)
  end
end

"""
second derivative of the drainage capillary pressure with respect to water saturation
dpc_drain(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc0::Real=1.0e6)
"""
function d2pc_drain3(sw::Real, pce::Real, swc::Real; labda::Real=2.4, pc_star::Real=1.0e6, c::Real=1e-3)
  sw_star=(1-swc)*exp(-labda*log(pc_star/pce))+swc
  b = -pce*c/(labda*(1-swc))*((sw_star-swc)/(1-swc))^(-1.0/labda-1.0)
  a = pc_star + b*exp(c)
  pc0 = a-b*exp(sw_star+c)
  if sw>sw_star
    res=-1.0/((1-swc)^2*labda)*(-1.0/labda-1.0)*pce*((sw-swc)/(1-swc))^(-1.0/labda-2.0)
  elseif 0.0<=sw<=sw_star
    res=b/(sw_star-sw+c)^2
  else
    res=b/(sw_star+c)^2
  end
end

"""
Imbibition capillary pressure curve
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e6, pc02::Real=1.0e6)
"""
function pc_imb3(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=1.0, pc_star1::Real=1.0e6, pc_star2::Real=1.0e6, c::Real=1e-3)
  pc1=pc_drain3(sw, pce, swc, labda=labda, pc_star=pc_star1, c=c)
  pc2=pc_drain3(1-sw, pce, sor, labda=labda, pc_star=pc_star2, c=c)
  return (0.5*(1+cos(teta)))^b*pc1-(0.5*(1-cos(teta)))^b*pc2
end

"""
derivative of the imbibition capillary pressure curve w.r.t water saturation
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e6, pc02::Real=1.0e6)
"""
function dpc_imb3(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda1::Real=2.4, labda2::Real=2.4, b::Real=1.0, pc_star1::Real=1.0e6, pc_star2::Real=1.0e6, c::Real=1e-3)
  dpc1=dpc_drain3(sw, pce, swc, labda=labda1, pc_star=pc_star1, c=c)
  dpc2=dpc_drain3(1-sw, pce, sor, labda=labda2, pc_star=pc_star2, c=c)
  return (0.5*(1+cos(teta)))^b*dpc1+(0.5*(1-cos(teta)))^b*dpc2
end

"""
derivative of the imbibition capillary pressure curve w.r.t water saturation
pc_imb(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda::Real=2.4, b::Real=0.6, pc01::Real=1.0e6, pc02::Real=1.0e6)
"""
function d2pc_imb3(sw::Real, pce::Real, swc::Real, sor::Real;
  teta::Real=0.785, labda1::Real=2.4, labda2::Real=2.4, b::Real=1.0, pc_star1::Real=1.0e6, pc_star2::Real=1.0e6, c::Real=1e-3)
  dpc1=d2pc_drain3(sw, pce, swc, labda=labda1, pc_star=pc_star1, c=c)
  dpc2=d2pc_drain3(1-sw, pce, sor, labda=labda2, pc_star=pc_star2, c=c)
  return (0.5*(1+cos(teta)))^b*dpc1+(0.5*(1-cos(teta)))^b*dpc2
end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng