text:
using Base: cumprod
using TensorCast
using NNlib

"""
    contentaddress(key::AbstractArray{T, 3}, mem::AbstractArray{T, 3}, β::AbstractArray{S, 2}, K=weightedcosinesim) where {T, S}

Compute the cosine similarity between all rows of memory M and the key k.
β acts as sharpener: high values concentrate weights.

# Arguments
- `k`: (W x R x B) keys
- `M`: (N x W x B) memory
- `β`: (R x B) strengths
"""

function contentaddress(k, M, β)
    dot = batched_mul(M, k)
    M = M .+ eps(eltype(M))
    k = k .+ eps(eltype(k))
    norm_k = sum(k.^2, dims=1)
    norm_M = sum(M.^2, dims=2)
    norm = sqrt.(batched_mul(norm_M, norm_k))
    β = reshape(β, 1, size(β)...)
    weightedsimilarity = dot .* β ./ norm
    softmax(weightedsimilarity; dims=1)
end

"""
    memoryretention(wr, f)

# Arguments
- `wr` - (N, R, B): previous iteration read weights.
- `f` - (R, B): free gates
Determine how much each memory location will not be freed by the free gates.
Returns a tensor of size (N x B)
"""
function memoryretention(wr, f)
    rs = one(eltype(wr)) .- wr.*reshape(f, 1, size(f)...)
    rs = prod(rs; dims=2)
    return dropdims(rs; dims=2)
end

_usage(u_prev, ww_prev, 𝜓) = (u_prev + ww_prev - (u_prev.*ww_prev)) .* 𝜓

"""
usage(u_prev, ww_prev, wr_prev, freegate)

# Arguments
- `u_prev` - (N, B) usage after previous iteration
- `ww_prev` - (N, 1, B) previous write weights
- `wr_prev` - (N, R, B) previous read weights
- `freegate` - (R, B)

Calculate the usage vector of the memory rows.
A row is considered used (``u[i]=1``) if they have recently been written to and haven't been retained by the free gates (``𝜓[i]=1``)
"""
function usage(u_prev, ww_prev, wr_prev, freegate)
    if ndims(ww_prev) == 3
        ww_prev = dropdims(ww_prev; dims=2)
    end
    𝜓 = memoryretention(wr_prev, freegate)
    _usage(u_prev, ww_prev, 𝜓)
end

const _EPSILON = 1f-6

"""
    allocationweighting(u::AbstractMatrix; eps::AbstractFloat=_EPSILON)
# Arguments
- `u`: (N x B) usage tensor

# Returns
- `a`: (N x B) tensor

"""
function  allocationweighting(u::T; eps::AbstractFloat=_EPSILON) where T<:AbstractMatrix
    u = eps .+ (1-eps) .* u
    ϕ = [sortperm(u[:, i]) for i in 1:size(u, 2)]
    ϕ = reshape(vcat(ϕ...), size(u))
    for i in CartesianIndices(ϕ)
        ϕ[i] += size(ϕ, 1)*(i[2]-1)
    end
    sortedusage = u[ϕ]
    prodsortedusage = cumprodexclusive(sortedusage;dims=1)
    sortedalloc = (1 .- sortedusage) .* prodsortedusage
    a = sortedalloc[ϕ]
    a::T
end

using Zygote: @adjoint
# The sorting of allocation weighting introduce discontinuities
# in the backward pass, so we set the pullback to 1
@adjoint allocationweighting(u::AbstractMatrix; eps=_EPSILON) =
allocationweighting(u; eps=eps), Δ -> (Δ, )

"""
    writeweight(cw::AbstractArray{T, 3},
                 a::AbstractArray{T, 2},
                gw::AbstractArray{T, 1},
                ga::AbstractArray{T, 1}) where T
Batch version
#Arguments
- `cw`: (N x 1 x B) write content weighting
- `a`: (N x B) allocation weighting
- `gw`: (1 x B) write gate
- `ga`: (1 x B) allocation gate
a:
"""
function writeweight(cw::AbstractArray{T, 3},
                     a::AbstractArray{T, 2},
                     gw::AbstractArray{T, 2},
                     ga::AbstractArray{T, 2}) where T
    a = reshape(a, size(a,1), 1, size(a, 2))
    gw = reshape(gw, 1, size(gw)...)
    ga = reshape(ga, 1, size(ga)...)
    return gw.*(ga .* a + (one(T) .- ga).*cw)
end

"""
    precedenceweight(p_prev::AbstractArray{T, 2}, ww::AbstractArray{T, 3}) where T

p_prev: (N, B)
ww: (N, 1, B)

return p; (N, B)
"""
function precedenceweight(p_prev::AbstractArray{T, 2}, ww::AbstractArray{T, 3}) where T
    ww = dropdims(ww; dims=2)
    wwsum = sum(ww; dims=1) #(1 x B)
    return (one(T).-wwsum).*p_prev + ww
end

"""
    updatelinkmatrix(L::AbstractArray{T, 3}, p::AbstractArray{T, 2}, ww::AbstractArray{T, 3}) where T

# Arguments
- `L`: (N x N x B)
- `p`: (N x B)
- `ww`: (N x 1 x B)
"""
function updatelinkmatrix!(L::AbstractArray{T, 3}, p::AbstractArray{T, 2}, ww::AbstractArray{T, 3}) where T
    ww = dropdims(ww; dims=2)
    N, _, B = size(L)
    for b in 1:B
        for i in 1:N
            for j in 1:N
                if i != j
                    L[i, j, b] = (one(T) - ww[i, b] - ww[j, b]) * L[i, j, b] + ww[i, b]*p[j, b]
                end
            end
        end
    end
    L
end

"""
    forwardweight(L::AbstractArray{T, 3}, wr::AbstractArray{T, 3}) where T
Location weight for reading the next-written-to location.

# Arguments
- `L`: (N x N x B) link matrix
- `wr`: (N x R x B) read weights

See also: [`backwardweight`](@ref)
"""
function forwardweight(L::AbstractArray{T, 3}, wr::AbstractArray{T, 3}) where T
    batched_mul(L, wr)
end

"""
    backwardweight(L::AbstractArray{T, 3}, wr::AbstractArray{T, 3}) where T
Location weight for reading the previous-written-to location.

# Arguments
- `L`: (N x N x B) link matrix
- `wr`: (N x R x B) read weights

See also: [`forwardweight`](@ref)
"""
function backwardweight(L::AbstractArray{T, 3}, wr::AbstractArray{T, 3}) where T
    batched_mul(PermutedDimsArray(L, (2, 1, 3)), wr)
end

"""
    readweigth(backw::AbstractArray::{T, 2},
    cr::AbstractArray{T, 3},
    forw::AbstractArray{T, 2},
    readmode::AbstractArray{T, 2}) where T

Interpolate the backward weighting, content weighting and forward weighting.
readmode is a vector of size 3 summing to 1.

# Arguments
- `backw`, `cr`, `forw`: (N x R x B)
- `readmode`: (3 x R x B)
```
# Returns
- (N x R x B) tensor represented each read heads readweights
"""
function readweight(backw::AbstractArray{T, 3}, cr::AbstractArray{T, 3}, forw::AbstractArray{T, 3}, readmode::AbstractArray{T, 3}) where T
    @cast out[n, r, b] := backw[n, r, b]*readmode[1, r, b] + cr[n, r, b]*readmode[2, r, b] + forw[n, r, b] * readmode[3, r, b]
end


text:
## Exercise 1-1
## Whenever you are experimenting with a new feature, you should try to make mistakes. For example, in the “Hello, World!” program, what happens if you leave out one of the quotation marks? What if you leave out both? What if you spell println wrong?

## This kind of experiment helps you remember what you read; it also helps when you are programming, because you get to know what the error messages mean. It is better to make mistakes now and on purpose rather than later and accidentally.

## 1. In a print statement, what happens if you leave out one of the parentheses, or both?
println("Ans 1:")
# println("wello"  # error, syntax incomplete
# println "wello"  # error syntax

## 2. If you are trying to print a string, what happens if you leave out one of the quotation marks, or both?
println("Ans 2:")
# println("wello)  # error, syntax incomplete
# println(wello)  # error, UndefVarError

## 3. You can use a minus sign to make a negative number like -2. What happens if you put a plus sign before a number? What about 2++2?
println("Ans 3:")
println(-2) # works
println(2 + +2) # works
# println(2++2)  # Error, UndefVarError

## 4. In math notation, leading zeros are okay, as in 02. What happens if you try this in Julia?
println("Ans 4:")
println(02) # works
println(02.)  # 2.0
println(2.0000000)

## 5. What happens if you have two values with no operator between them?
println("Ans 5:")
# println(3 5)  # error, syntax, missing , or ) in argument list

println("End.")

ext: jl


text:
abstract type AbstractPowerModel end

"""
    PowerModelConstantCp(cp)

Models will assume a constant cp value as provided

# Arguments
- `cp::Float`: constant power coefficient value
- 'pp::TI': exponent for adjusting power for wind turbine yaw
"""
struct PowerModelConstantCp{TF,TI} <: AbstractPowerModel
    cp::TF
    pp::TI
end
PowerModelConstantCp(x) = PowerModelConstantCp(x,2)

"""
    PowerModelCpPoints(vel_points, cp_points)

Models will use adjust cp based on cp curve using linear interpolation of provided points

# Arguments
- `vel_points::Array{N,Float}`: wind speed values in m/s
- `cp_points::Array{N,Float}`: power coefficient values corresponding to the provided speeds
- 'pp::TF': exponent for adjusting power for wind turbine yaw
"""
struct PowerModelCpPoints{ATF,TF} <: AbstractPowerModel
    vel_points::ATF
    cp_points::ATF
    pp::TF
end
PowerModelCpPoints(x,y) = PowerModelCpPoints(x, y, 2)

"""
    PowerModelPowerPoints(vel_points, cp_points)

Models will use adjust wind turbine power based on power curve using linear interpolation of
provided points

# Arguments
- `vel_points::Array{N,Float}`: wind speed values in m/s
- `power_points::Array{N,Float}`: power values corresponding to the provided speeds
- 'pp::TF': exponent for adjusting power for wind turbine yaw
"""
struct PowerModelPowerPoints{ATF, TF} <: AbstractPowerModel
    vel_points::ATF
    power_points::ATF
    pp::TF
end
PowerModelPowerPoints(x,y) = PowerModelPowerPoints(x, y, 2)

"""
    PowerModelPowerCurveCubic()

Power will be calculated based on turbine specifications assuming a cubic power curve. Note
that this method is inherently incorrect and should only be used for theoretical purposes
or after careful validation.

# Arguments
- 'pp::TF': exponent for adjusting power for wind turbine yaw

"""
struct PowerModelPowerCurveCubic{TF} <: AbstractPowerModel
    pp::TF
end
PowerModelPowerCurveCubic() = PowerModelPowerCurveCubic(2)
"""
    calculate_power_from_cp(generator_efficiency, air_density, rotor_area, cp, wt_velocity)

Calculate the power for a wind turbine based on standard theory for region 2

# Arguments
- `generator_efficiency::Float`: Efficiency of the turbine generator
- `air_density::Float`: Air density
- `rotor_area::Float`: Rotor-swept area of the wind turbine
- `cp::Float`: Power coefficient of the wind turbine
- `wt_velocity::Float`: Inflow velocity to the wind turbine
"""
function calculate_power_from_cp(generator_efficiency, air_density, rotor_area, cp, wt_velocity, wt_yaw; pp=2)
    power = generator_efficiency*(0.5*air_density*rotor_area*cp*(cos(wt_yaw)^pp)*wt_velocity^3)
    return power
end

"""
    calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, power_model)

Calculate the power for a wind turbine based on standard theory for region 2 using a constant cp

# Arguments
- `generator_efficiency::Float`: Efficiency of the turbine generator
- `air_density::Float`: Air density
- `rotor_area::Float`: Rotor-swept area of the wind turbine
- `wt_velocity::Float`: Inflow velocity to the wind turbine
- `cut_in_speed::Float`: cut in speed of the wind turbine
- `rated_speed::Float`: rated speed of the wind turbine
- `cut_out_speed::Float`: cut out speed of the wind turbine
- `rated_power::Float`: rated power of the wind turbine
- `power_model::PowerModelConstantCp`: Struct containing the cp value to be used in region 2
"""
function calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, wt_yaw, cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model::PowerModelConstantCp)

    if wt_velocity < cut_in_speed
        power = 0.0
    elseif wt_velocity < rated_speed
        # extract cp_value
        cp = power_model.cp
        pp = power_model.pp
        power = calculate_power_from_cp(generator_efficiency, air_density, rotor_area, cp, wt_velocity, wt_yaw, pp=pp)
        if power > rated_power
            power = rated_power
        end
    elseif wt_velocity < cut_out_speed
        power = rated_power
    elseif wt_velocity > cut_out_speed
        power = 0.0
    end

    return power

end

"""
    calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model)

Calculate the power for a wind turbine based on a cp curve with linear interpolation

# Arguments
- `generator_efficiency::Float`: Efficiency of the turbine generator
- `air_density::Float`: Air density
- `rotor_area::Float`: Rotor-swept area of the wind turbine
- `wt_velocity::Float`: Inflow velocity to the wind turbine
- `cut_in_speed::Float`: cut in speed of the wind turbine
- `rated_speed::Float`: rated speed of the wind turbine
- `cut_out_speed::Float`: cut out speed of the wind turbine
- `rated_power::Float`: rated power of the wind turbine
- `power_model::PowerModelCpPoints`: Struct containing the velocity and cp values defining the cp curve
"""
function calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, wt_yaw,
    cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model::PowerModelCpPoints)

    # obtain exponent for adjusting cp to yaw
    pp = power_model.pp

    # initialize power to zero
    power = 0.0

    # use specs if inflow wind speed is less than the wind speeds provided in the power curve
    if wt_velocity < power_model.vel_points[1]

        # calculated wind turbine power
        if wt_velocity < cut_in_speed
            power = 0.0
        elseif wt_velocity < rated_speed
            # use cp value corresponding to lowest provided velocity point
            cp = power_model.cp_points[1]

            # calculate power
            power = calculate_power_from_cp(generator_efficiency, air_density, rotor_area, cp, wt_velocity, wt_yaw, pp=pp)
        end

    # use cp points where provided
    elseif wt_velocity < power_model.vel_points[end]

        # estimate cp_value using linear interpolation
        cp = linear(power_model.vel_points, power_model.cp_points, wt_velocity)

        # calculate power
        power = calculate_power_from_cp(generator_efficiency, air_density, rotor_area, cp, wt_velocity, wt_yaw, pp=pp)

    # use specs if above vel_points max
    else

        if wt_velocity <= cut_out_speed
            # use cp value corresponding to highest provided velocity point
            cp = power_model.cp_points[end]

            # calculate power
            power = calculate_power_from_cp(generator_efficiency, air_density, rotor_area, cp, wt_velocity, wt_yaw, pp=pp)
        end

    end

    return power

end

"""
    calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model)

Calculate the power for a wind turbine based on a pre-determined power curve with linear
    interpolation

# Arguments
- `generator_efficiency::Float`: Efficiency of the turbine generator
- `air_density::Float`: Air density
- `rotor_area::Float`: Rotor-swept area of the wind turbine
- `wt_velocity::Float`: Inflow velocity to the wind turbine
- `cut_in_speed::Float`: cut in speed of the wind turbine
- `rated_speed::Float`: rated speed of the wind turbine
- `cut_out_speed::Float`: cut out speed of the wind turbine
- `rated_power::Float`: rated power of the wind turbine
- `power_model::PowerModelPowerPoints`: Struct containing the velocity and power values
    defining the power curve
"""
function calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, wt_yaw,
    cut_in_speed, rated_speed, cut_out_speed, rated_power,
    power_model::PowerModelPowerPoints; return_derivative=false)

    # get exponent for yaw adjustment
    pp = power_model.pp

    power = 0.0
    dpower_dvelocity = 0.0
    # use specs if inflow wind speed is less than the wind speeds provided in the power curve
    if wt_velocity < power_model.vel_points[1]

        # calculated wind turbine power
        if wt_velocity < rated_speed
            # use power value corresponding to lowest provided velocity point
            power = (cos(wt_yaw)^pp)*linear([cut_in_speed, power_model.vel_points[1]], [0.0, power_model.power_points[1]], wt_velocity)
        end

    # use power points where provided
    elseif wt_velocity < power_model.vel_points[end]

        # calculate power
        power = (cos(wt_yaw)^pp)*linear(power_model.vel_points, power_model.power_points, wt_velocity)
        if return_derivative
            dpower_dvelocity = (cos(wt_yaw)^pp)*gradient(power_model.vel_points, power_model.power_points, wt_velocity)
        end
    # use specs if above vel_points max
    else

        if wt_velocity <= cut_out_speed
            # use power corresponding to highest wind speed provided
            power = (cos(wt_yaw)^pp)*power_model.power_points[end]
        end
    end

    if !return_derivative
        return power
    else
        return power, dpower_dvelocity
    end
end

"""
    calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model)

Calculates wind turbine power using a cubic estimation based on turbine specifications
    as defined in https://github.com/byuflowlab/iea37-wflo-casestudies/blob/master/cs3-4/iea37-cs3-announcement.pdf

# Arguments
- `generator_efficiency::Float`: Efficiency of the turbine generator
- `air_density::Float`: Air density
- `rotor_area::Float`: Rotor-swept area of the wind turbine
- `wt_velocity::Float`: Inflow velocity to the wind turbine
- `cut_in_speed::Float`: cut in speed of the wind turbine
- `rated_speed::Float`: rated speed of the wind turbine
- `cut_out_speed::Float`: cut out speed of the wind turbine
- `rated_power::Float`: rated power of the wind turbine
- `power_model::PowerModelPowerCurveCubic`: Empty struct
"""
function calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, wt_yaw,
    cut_in_speed, rated_speed, cut_out_speed, rated_power,
    power_model::PowerModelPowerCurveCubic)

    pp = power_model.pp

    if wt_velocity < cut_in_speed
        power = 0.0
    elseif wt_velocity < rated_speed
        power = rated_power*((wt_velocity - cut_in_speed)/(rated_speed - cut_in_speed))^3
        # power = rated_power*((wt_velocity)/(rated_speed))^3
    elseif wt_velocity < cut_out_speed
        power = rated_power
    elseif wt_velocity > cut_out_speed
        power = 0.0
    end

    return power*(cos(wt_yaw)^pp)

end

"""
    calculate_turbine_power(generator_efficiency, cut_in_speed, cut_out_speed, rated_speed,
    rated_power, rotor_diameter, wt_velocity, power_model::AbstractPowerModel, air_density)

Calculate the power for all wind turbines. Dispaches to desired power model.

# Arguments
- `generator_efficiency::Array{Float,nTurbines}`
- `cut_in_speed::Array{Float,nTurbines}`
- `cut_out_speed::Array{Float,nTurbines}`
- `rated_speed::Array{Float,nTurbines}`
- `rated_power::Array{Float,nTurbines}`
- `rotor_diameter::Array{Float,nTurbines}`
- `wt_velocity::Array{Float,nTurbines}`: turbine effective wind speeds for current state only
- `power_model::AbstractPowerModel)
- `air_density::Float`
"""
function calculate_turbine_power(generator_efficiency, cut_in_speed, cut_out_speed,
    rated_speed, rated_power, rotor_diameter, wt_velocity, wt_yaw, power_model::AbstractPowerModel,
    air_density; return_derivatives=false)

    # calculated wind turbine rotor-swept area
    rotor_area = pi*(rotor_diameter^2)/4.0

    if !return_derivatives
        wt_power = calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, wt_yaw, cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model)
        return wt_power
    else
        wt_power, dp_du = calculate_power(generator_efficiency, air_density, rotor_area, wt_velocity, wt_yaw, cut_in_speed, rated_speed, cut_out_speed, rated_power, power_model, return_derivative=return_derivatives)
        return wt_power, dp_du
    end
end

"""
    turbine_powers_one_direction((generator_efficiency, cut_in_speed, cut_out_speed,
        rated_speed, rated_power, rotor_diameter, turbine_inflow_velcities, air_density, power_model::AbstractPowerModel)

Calculate the power for all wind turbines for a given state

# Arguments
- `generator_efficiency::Array{Float,nTurbines}`
- `cut_in_speed::Array{Float,nTurbines}`
- `cut_out_speed::Array{Float,nTurbines}`
- `rated_speed::Array{Float,nTurbines}`
- `rated_power::Array{Float,nTurbines}`
- `rotor_diameter::Array{Float,nTurbines}`
- `turbine_inflow_velcities::Array{Float,nTurbines}`: for current state only
- `air_density::Float`
- `power_models::Array{nturbines})` elements of array should be be of sub-types or AbstractPowerModel
"""
function turbine_powers_one_direction(generator_efficiency, cut_in_speed, cut_out_speed,
    rated_speed, rated_power, rotor_diameter, turbine_inflow_velcities, turbine_yaw, air_density,
    power_models; jac=nothing)

    # get number of turbines and rotor sample point
    nturbines = length(rotor_diameter)

    arr_type = promote_type(typeof(generator_efficiency[1]),typeof(cut_in_speed[1]),typeof(cut_out_speed[1]),typeof(rated_speed[1]),
                            typeof(rated_power[1]),typeof(rotor_diameter[1]),typeof(turbine_inflow_velcities[1]),typeof(turbine_yaw[1]))
    wt_power = zeros(arr_type, nturbines)

    if jac === nothing
        for d=1:nturbines
            wt_power[d] = calculate_turbine_power(generator_efficiency[d], cut_in_speed[d], cut_out_speed[d], rated_speed[d], rated_power[d], rotor_diameter[d], turbine_inflow_velcities[d], turbine_yaw[d], power_models[d], air_density)
        end
    else
        for d=1:nturbines
            wt_power[d], jac[d] = calculate_turbine_power(generator_efficiency[d], cut_in_speed[d], cut_out_speed[d], rated_speed[d], rated_power[d], rotor_diameter[d], turbine_inflow_velcities[d], turbine_yaw[d], power_models[d], air_density; return_derivatives=true)
        end
    end
    return wt_power
end

"""
    calculate_state_aeps(turbine_x, turbine_y, turbine_z, rotor_diameter,
    hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
    cut_out_speed, rated_speed, rated_power, wind_resource, power_models::Array{AbstractPowerModel}, model_set::AbstractModelSet;
    rotor_sample_points_y=[0.0], rotor_sample_points_z=[0.0])

Calculate AEP for each requested state respectively

# Arguments
- `turbine_x::Array{Float,nTurbines}`: turbine east-west locations in the global
    reference frame
- `turbine_y::Array{Float,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_z::Array{Float,nTurbines}`: turbine base height in the global reference frame
- `rotor_diameter::Array{Float,nTurbines}`
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `ct_model::AbstractThrustCoefficientModel`: defines how the thrust coefficient changes
    with state etc
- `generator_efficiency::Array{Float,nTurbines}`
- `cut_in_speed::Array{Float,nTurbines}`
- `cut_out_speed::Array{Float,nTurbines}`
- `rated_speed::Array{Float,nTurbines}`
- `rated_power::Array{Float,nTurbines}`
- `wind_resource::DiscretizedWindResource`: wind resource discreption (directions, speeds,
    frequencies, etc)
- `power_model::Array{nTurbines}`: elemenst of array should be sub-types of AbstractPowerModel
- `model_set::AbstractModelSet`: defines wake-realated models to be used in analysis
- rotor_sample_points_y::Array{TF,N}`: horizontal wind location of points to sample across
    the rotor swept area when calculating the effective wind speed for the wind turbine.
    Points are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- rotor_sample_points_z::Array{TF,N}`: vertical wind location of points to sample across the
    rotor swept area when calculating the effective wind speed for the wind turbine. Points
    are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- `hours_per_year::Float`: hours per year (averaged for leap year by default)
"""
function calculate_state_aeps(turbine_x, turbine_y, turbine_z, rotor_diameter,
            hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
            cut_out_speed, rated_speed, rated_power, wind_resource, power_models, model_set::AbstractModelSet;
            rotor_sample_points_y=[0.0], rotor_sample_points_z=[0.0], hours_per_year=365.25*24.0)

    wind_probabilities = wind_resource.wind_probabilities

    nstates = length(wind_probabilities)

    arr_type = promote_type(typeof(turbine_x[1]),typeof(turbine_y[1]),typeof(turbine_z[1]),typeof(rotor_diameter[1]),typeof(hub_height[1]),typeof(turbine_yaw[1]),
                typeof(generator_efficiency[1]),typeof(cut_in_speed[1]),typeof(cut_out_speed[1]),typeof(rated_speed[1]),typeof(rated_power[1]))
    state_energy = zeros(arr_type,nstates)

    for i = 1:nstates

        rot_x, rot_y = rotate_to_wind_direction(turbine_x, turbine_y, wind_resource.wind_directions[i])

        sorted_turbine_index = sortperm(rot_x)

        turbine_velocities = turbine_velocities_one_direction(rot_x, rot_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
                            sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                            model_set, wind_farm_state_id=i, velocity_only=true)

        wt_power = turbine_powers_one_direction(generator_efficiency, cut_in_speed, cut_out_speed, rated_speed,
                            rated_power, rotor_diameter, turbine_velocities, turbine_yaw, wind_resource.air_density, power_models)

        state_power = sum(wt_power)
        state_energy[i] = state_power*hours_per_year*wind_probabilities[i]
    end

    return state_energy
end

function calculate_state_aep(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height,
    turbine_yaw, ct_model, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed,
    rated_power, power_models, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
    model_set; wind_farm_state_id=1, hours_per_year=365.25*24.0, velocity_only=true)

    rot_x, rot_y = rotate_to_wind_direction(turbine_x, turbine_y, wind_resource.wind_directions[wind_farm_state_id])

    sorted_turbine_index = sortperm(rot_x)

    turbine_velocities = turbine_velocities_one_direction(rot_x, rot_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
                        sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                        model_set, wind_farm_state_id=wind_farm_state_id, velocity_only=true)

    wt_power = turbine_powers_one_direction(generator_efficiency, cut_in_speed, cut_out_speed, rated_speed,
                        rated_power, rotor_diameter, turbine_velocities, turbine_yaw, wind_resource.air_density, power_models)

    state_power = sum(wt_power)

    state_energy = state_power*hours_per_year*wind_resource.wind_probabilities[wind_farm_state_id]

    return state_energy

end

"""
    calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
    hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
    cut_out_speed, rated_speed, rated_power, wind_resource, power_models::Array{AbstractPowerModel}, model_set::AbstractModelSet;
    rotor_sample_points_y=[0.0], rotor_sample_points_z=[0.0], hours_per_year=365.25*24.0)

Calculate wind farm AEP

# Arguments
- `turbine_x::Array{Float,nTurbines}`: turbine east-west locations in the global
    reference frame
- `turbine_y::Array{Float,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_z::Array{Float,nTurbines}`: turbine base height in the global reference frame
- `rotor_diameter::Array{Float,nTurbines}`
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `ct_model::AbstractThrustCoefficientModel`: defines how the thrust coefficient changes
    with state etc
- `generator_efficiency::Array{Float,nTurbines}`
- `cut_in_speed::Array{Float,nTurbines}`
- `cut_out_speed::Array{Float,nTurbines}`
- `rated_speed::Array{Float,nTurbines}`
- `rated_power::Array{Float,nTurbines}`
- `wind_resource::DiscretizedWindResource`: wind resource discreption (directions, speeds,
    frequencies, etc)
- `power_model::Array{)`: elements of array should be sub types of AbstractPowerModel
- `model_set::AbstractModelSet`: defines wake-realated models to be used in analysis
- `rotor_sample_points_y::Array{TF,N}`: horizontal wind location of points to sample across
    the rotor swept area when calculating the effective wind speed for the wind turbine.
    Points are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- `rotor_sample_points_z::Array{TF,N}`: vertical wind location of points to sample across the
    rotor swept area when calculating the effective wind speed for the wind turbine. Points
    are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- `hours_per_year::Float`: hours per year (averaged for leap year by default)
"""
function calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
            hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
            cut_out_speed, rated_speed, rated_power, wind_resource, power_models, model_set::AbstractModelSet;
            rotor_sample_points_y=[0.0], rotor_sample_points_z=[0.0], hours_per_year=365.25*24.0, distributed=false)

    # find how many wind states are being calculated
    nstates = length(wind_resource.wind_probabilities)

    # state_energy = Vector{typeof(wind_farm.turbine_x[1])}(undef,nstates)
    arr_type = promote_type(typeof(turbine_x[1]),typeof(turbine_y[1]),typeof(turbine_z[1]),typeof(rotor_diameter[1]),typeof(hub_height[1]),typeof(turbine_yaw[1]),
                typeof(generator_efficiency[1]),typeof(cut_in_speed[1]),typeof(cut_out_speed[1]),typeof(rated_speed[1]),typeof(rated_power[1]))

    # calculate AEP in parallel using multi-threading
    if Threads.nthreads() > 1
        state_aep = zeros(arr_type,nstates)
        Threads.@threads for i = 1:nstates

            state_aep[i] = calculate_state_aep(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height,
                turbine_yaw, ct_model, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed,
                rated_power, power_models, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                model_set; wind_farm_state_id=i, hours_per_year=hours_per_year, velocity_only=true)

        end

        AEP = sum(state_aep)

    # calculate AEP in serial or in parallel using distributed processing
    else
        AEP = @sync @distributed (+) for i = 1:nstates
            state_aep = calculate_state_aep(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height,
                turbine_yaw, ct_model, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed,
                rated_power, power_models, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                model_set; wind_farm_state_id=i, hours_per_year=hours_per_year, velocity_only=true)
        end
    end

    return AEP
end


text:
abstract type EngineMode end

mutable struct PhysicsState{T<:Union{MJSim }}
    model::T
    pert::RefValue{mjvPerturb}
    elapsedsim::Float64
    timer::RateTimer
    lock::ReentrantLock

    function PhysicsState(model::Union{MJSim })
        pert = Ref(mjvPerturb())
        mjv_defaultPerturb(pert)
        new{typeof(model)}(model, pert, 0, RateTimer(), ReentrantLock())
    end
end

Base.@kwdef mutable struct UIState
    scn::RefValue{mjvScene} = Ref(mjvScene())
    cam::RefValue{mjvCamera} = Ref(mjvCamera())
    vopt::RefValue{mjvOption} = Ref(mjvOption())
    con::RefValue{mjrContext} = Ref(mjrContext())
    figsensor::RefValue{mjvFigure} = Ref(mjvFigure())

    showinfo::Bool = true
    showsensor::Bool = false
    speedmode::Bool = false
    speedfactor::Float64 = 1 / 10

    reversed::Bool = false
    paused::Bool = true
    shouldexit::Bool = false

    reward::Float64 = 0

    lastrender::Float64 = 0
    refreshrate::Float64 = 0
    realtimerate::Float64 = 0
    io1::IOBuffer = IOBuffer()
    io2::IOBuffer = IOBuffer()

    lock::ReentrantLock = ReentrantLock()
end

mutable struct Engine{T,M}
    phys::PhysicsState{T}
    ui::UIState
    mngr::WindowManager
    handlers::Vector{EventHandler}

    modes::M
    modehandlers::Vector{EventHandler}
    curmodeidx::Int

    ffmpeghandle::Maybe{Base.Process}
    framebuf::Vector{UInt8}
    videodst::Maybe{String}
    min_refreshrate::Int

    function Engine(windowsize::NTuple{2,Integer}, model::Union{MJSim }, modes::Tuple{Vararg{EngineMode}})
        window = create_window(windowsize..., "MUJOCO_SIM")
        try
            phys = PhysicsState(model)
            ui = UIState()
            mngr = WindowManager(window)

            mjv_defaultScene(ui.scn)
            mjv_defaultCamera(ui.cam)
            mjv_defaultOption(ui.vopt)
            mjr_defaultContext(ui.con)
            mjv_defaultFigure(ui.figsensor)

            sim = getsim(model)
            mjv_makeScene(sim.m, ui.scn, MAXGEOM) # TODO calculate MAXGEOM form model and moes
            mjr_makeContext(sim.m, ui.con, FONTSCALE)

            alignscale!(ui, sim)
            init_figsensor!(ui.figsensor)

            e = new{typeof(model),typeof(modes)}(
                phys,
                ui,
                mngr,
                EventHandler[],

                modes,
                handlers(ui, phys, first(modes)),
                1,

                nothing,
                UInt8[],
                nothing,
                min(map(GetRefreshRate, GLFW.GetMonitors())..., MIN_REFRESHRATE),
            )

            e.handlers = handlers(e)
            register!(mngr, e.handlers...)

            return e
        catch e
            GLFW.DestroyWindow(window)
            rethrow(e)
        end
    end
end


text:
using Distributions
import StatsBase.ecdf
#using Plots

function s_weibull( x::Float64, k::Float64=1.0, λ::Float64=1.0)
    return 1.0-cdf( Weibull(k, λ), x)
end

#using PyPlot
# figure()
# plot( collect(0.:1.0/100.:5.),[s_weibull(i, 2.2, 1.75) for i in collect(0.:1.0/100.:5.)], linewidth=1.1, label = "\$ Y\\^{(0)} survival with \$l=1\$"  )
# plot( collect(0.:1.0/100.:5.),[s_weibull(i, 1.2, 2.2) for i in collect(0.:1.0/100.:5.)], linewidth=1.1, label = "\$ Y^{(1)} survival with \$l=0\$"   )
# plot( collect(0.:1.0/100.:5.),[s_weibull(i, 5.3, 1.5) for i in collect(0.:1.0/100.:5.)], linewidth=1.1, label = "\$ Y^{(2)} survival with \$l=0\$" )
# ylabel("\$ S(t) \$")
# xlabel("\$ t \$")
# title("Simulated data true survival functions")
# legend(loc=1)
# savefig("weibulls.png")

c_sw = [5.0,2.2]
c_ws = [1.5,3.4]
c_m = [2.1,2.0]

function s_comprisk_weibull( x::Float64, k₁::Float64=1.0, k₂::Float64=1.0,
                        λ₁::Float64=1.0,  λ₂::Float64=1.0)
    return exp( log( ccdf( Weibull(k₁, λ₁), x) ) + log( ccdf( Weibull(k₂, λ₂), x) ) )
end

# plot( collect(0.:1.0/100.:5.),[s_comprisk_weibull(i, 5.0, 1.2, 1.8, 2.2) for i in collect(0.:1.0/100.:5.)] )
# plot!( collect(0.:1.0/100.:5.),[s_comprisk_weibull(i, 5.0, 4.3, 1.8, 1.5) for i in collect(0.:1.0/100.:5.)] )

function sim_study_sampler_comprisk( n::Int64, δ::Float64, thr::Float64,
        k₀::Float64=1.0, λ₀::Float64=1.0, k₁::Float64=1.0, λ₁::Float64=1.0,
        k₂::Float64=1.0, λ₂::Float64=1.0,
        μ_co::Float64=1.0, σ_co::Float64=0.75, d_co::Float64=0.0, u_co::Float64=2.0)
    samp_Z = zeros( Float64, n)
    samp_δ = zeros( Float64, n)
    samp_Y = zeros( Float64, n)
    for i = 1:n
        samp_Z[i] = rand( TruncatedNormal( μ_co, σ_co, d_co, u_co) )
        if samp_Z[i] < thr
            samp_Y[i] = min( rand( Weibull( k₁, λ₁/(1.0-δ)^(1.0/k₁) ) ),
                                 rand( Weibull( k₀, λ₀/(δ)^(1.0/k₀) ) ) )
        else
            samp_Y[i] = min( rand( Weibull( k₂, λ₂/(1.0-δ)^(1.0/k₂) ) ),
                                 rand( Weibull( k₀, λ₀/(δ)^(1.0/k₀) ) ) )
        end
    end
    return [ samp_Y, samp_Z ]
end

function true_surv_comprisk( x::Float64, z::Float64, δ::Float64, thr::Float64,
    k₀::Float64=1.0, λ₀::Float64=1.0, k₁::Float64=1.0, λ₁::Float64=1.0,
    k₂::Float64=1.0, λ₂::Float64=1.0)
    if z <= thr
        return s_comprisk_weibull( x, k₁, k₀,
        λ₁/(1.0-δ)^(1.0/k₁), λ₀/(δ)^(1.0/k₀) )
    else
        return s_comprisk_weibull( x, k₂, k₀,
        λ₂/(1.0-δ)^(1.0/k₂), λ₀/(δ)^(1.0/k₀) )
    end
end

# Kaplan-Meier estiamtor function
function KMest(tiempo, muerte)
    t_KMest = sort( unique( tiempo ) )
    if t_KMest[1] != 0
        t_KMest = [0; t_KMest]
    end
    n_KMest = [count( i->(i >= j), tiempo) for j in t_KMest]
    d_KMest = [sum( muerte[ findall(j->j==i,tiempo) ] ) for i in t_KMest]
    s_KMest = cumprod( 1 .- d_KMest./n_KMest )
    return t_KMest, s_KMest
end

# plot( collect(0.:1.0/100.:5.), [true_surv_comprisk( i, 0.5, 0.1, 1.0, 5.0, 1.8, 1.2, 2.2, 4.3, 1.5) for i in collect(0.:1.0/100.:5.)])
# plot!( collect(0.:1.0/100.:5.), [true_surv_comprisk( i, 1.5, 0.1, 1.0, 5.0, 1.8, 1.2, 2.2, 4.3, 1.5) for i in collect(0.:1.0/100.:5.)])


text:
module Rotate

export @rotate, @rotaten

"""
    @rotate(x_1, x_2, ... x_n)

Rotate the values referred to by `x_i`.

Assign `x_2` to `x_1`, `x_3` to `x_2`, ..., and `x_1` to `x_n`.
Each argument is interpreted as the left hand side of an assignment.

`@rotate` is equivalent to `(x_1, x_2, ..., x_n) = (x_n, x_1, ..., x_{n-1})`

Compare to `circshift`.

## Example
```julia
@rotate(x, y, a[2])
```
"""
macro rotate(args...)
    return _rotate(1, args...)
end

"""
    @rotaten(n, x_1, x_2, ... x_n)

Like `@rotate`, but rotate by `n`, where `n` is an integer.
Use `n = -1` for circular shift left.
"""
macro rotaten(n, args...)
    return _rotate(n, args...)
end

function _rotate(n, args...)
    args_rhs = (circshift([args...], n)...,)
    expr = quote
        ($(esc.(args)...),) = ($(esc.(args_rhs)...),)
    end
    return expr
end

end # module


text:
module STREAMBenchmark

using CpuId, BenchmarkTools
using Statistics, Downloads
using Base.Threads: nthreads, @threads
using LoopVectorization

include("kernels.jl")
include("benchmarks.jl")
include("original.jl")

export memory_bandwidth, benchmark

end


text:
@naiveVariationalRule(:node_type     => AutoRegressiveMovingAverage,
                      :outbound_type => Message{GaussianMeanPrecision},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalARMAOutNPPPP)

@naiveVariationalRule(:node_type     => AutoRegressiveMovingAverage,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalARMAIn1PNPPP)

@naiveVariationalRule(:node_type     => AutoRegressiveMovingAverage,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalARMAIn2PPNPP)

@naiveVariationalRule(:node_type     => AutoRegressiveMovingAverage,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VariationalARMAIn3PPPNP)

@naiveVariationalRule(:node_type     => AutoRegressiveMovingAverage,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VariationalARMAIn4PPPPN)

max_stars_repo_name: wmkouw/ARMA

max_issues_repo_name: wmkouw/ARMA

max_forks_repo_name: wmkouw/ARMA


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.database;
#
# import javax.annotation.CheckReturnValue;
# import java.util.List;
#
# public interface DatabaseManager {
#
#     @CheckReturnValue
#     Database get(String name);
#
#     @CheckReturnValue
#     boolean contains(String name);
#     // TODO: Return type should be 'Database' but right now that would require 2 server calls in Cluster
#
#     void create(String name);
#
#     @CheckReturnValue
#     List<? extends Database> all();
#
#     interface Cluster extends DatabaseManager {
#
#         @Override
#         @CheckReturnValue
#         Database.Cluster get(String name);
#
#         @Override
#         @CheckReturnValue
#         List<Database.Cluster> all();
#     }
# }


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

@doc raw"""A label selector is a label query over a set of resources. The result of matchLabels and matchExpressions are ANDed. An empty label selector matches all objects. A null label selector matches no objects.

    IoK8sApimachineryPkgApisMetaV1LabelSelector(;
        matchExpressions=nothing,
        matchLabels=nothing,
    )

    - matchExpressions::Vector{IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement} : matchExpressions is a list of label selector requirements. The requirements are ANDed.
    - matchLabels::Dict{String, String} : matchLabels is a map of {key,value} pairs. A single {key,value} in the matchLabels map is equivalent to an element of matchExpressions, whose key field is \"key\", the operator is \"In\", and the values array contains only \"value\". The requirements are ANDed.
"""
mutable struct IoK8sApimachineryPkgApisMetaV1LabelSelector <: SwaggerModel
    matchExpressions::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement} } # spec name: matchExpressions
    matchLabels::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: matchLabels

    function IoK8sApimachineryPkgApisMetaV1LabelSelector(;matchExpressions=nothing, matchLabels=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1LabelSelector, Symbol("matchExpressions"), matchExpressions)
        setfield!(o, Symbol("matchExpressions"), matchExpressions)
        validate_property(IoK8sApimachineryPkgApisMetaV1LabelSelector, Symbol("matchLabels"), matchLabels)
        setfield!(o, Symbol("matchLabels"), matchLabels)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1LabelSelector

const _property_map_IoK8sApimachineryPkgApisMetaV1LabelSelector = Dict{Symbol,Symbol}(Symbol("matchExpressions")=>Symbol("matchExpressions"), Symbol("matchLabels")=>Symbol("matchLabels"))
const _property_types_IoK8sApimachineryPkgApisMetaV1LabelSelector = Dict{Symbol,String}(Symbol("matchExpressions")=>"Vector{IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement}", Symbol("matchLabels")=>"Dict{String, String}")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelector }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1LabelSelector))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelector }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1LabelSelector[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelector }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1LabelSelector[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1LabelSelector)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelector }, name::Symbol, val)
end


text:
using GaussianMixtureAlignment
using CoordinateTransformations
# Note: ProfileView and BenchmarkTools are not added by this package
using ProfileView
using BenchmarkTools

# small problem
xpts = [[0.,0.,0.], [3.,0.,0.,], [0.,4.,0.]]
ypts = [[1.,1.,1.], [1.,-2.,1.], [1.,1.,-3.]]
gmmx = IsotropicGMM([IsotropicGaussian(x, 1, 1) for x in xpts])
gmmy = IsotropicGMM([IsotropicGaussian(y, 1, 1) for y in ypts])

@btime gogma_align(gmmx, gmmy; maxblocks=1e3)
@btime tiv_gogma_align(gmmx, gmmy; maxblocks=1e3)

@ProfileView.profview gogma_align(gmmx, gmmy; maxblocks=1e5)
@ProfileView.profview tiv_gogma_align(gmmx, gmmy)

# steroid-sized problem
randpts = 10*rand(3,50)
randtform = AffineMap(10*rand(6)...)
gmmx = IsotropicGMM([IsotropicGaussian(randpts[:,i],1,1) for i=1:size(randpts,2)])
gmmy = randtform(gmmx)
min_overlap_score = overlap(gmmx,gmmx)
@btime res = gogma_align(gmmx,gmmy; maxblocks=1E3)
@btime tiv_res = tiv_gogma_align(gmmx,gmmy,0.5,0.5; maxblocks=1E3)

@ProfileView.profview res = gogma_align(gmmx, gmmy; maxblocks=1e4)
@ProfileView.profview tiv_res = tiv_gogma_align(gmmx, gmmy, 0.5, 0.5; maxblocks=1e4)


text:
using StatsModelComparisons, StanSample
using CSV, Printf, StatsPlots

ProjDir = @__DIR__

df = CSV.read(joinpath(ProjDir, "roachdata.csv"), DataFrame)
df.roach1 = df.roach1 / 100

roaches1_stan = "
data {
  int<lower=0> N;
  int<lower=0> K;
  vector[N] exposure2;
  vector[N] roach1;
  vector[N] senior;
  vector[N] treatment;
  int y[N];
}
transformed data {
  vector[N] log_expo;
  log_expo = log(exposure2);
}
parameters {
  vector[K] beta;
}
transformed parameters {
   vector[N] eta;
   eta = log_expo + beta[1] + beta[2] * roach1 + beta[3] * treatment
                + beta[4] * senior;
}
model {
  y ~ poisson_log(eta);
}
generated quantities {
  vector[N] log_lik;
  for (i in 1:N)
    log_lik[i] = poisson_log_lpmf(y[i] | eta[i]);
}
";

data = (N = size(df, 1), K = 4, y = Int.(df.y), roach1=df.roach1,
  exposure2=df.exposure2, senior=df.senior, treatment=df.treatment)
n = size(df, 1)

tmpdir = mktempdir()
sm1 = SampleModel("roaches1", roaches1_stan; tmpdir)
rc1 = stan_sample(sm1; data)

if success(rc1)
  stan_summary(sm1, true)
  nt1 = read_samples(sm1)

  # Compute LOO and standard error
  log_lik1 = nt1.log_lik'
  loo1, loos1, pk1 = psisloo(log_lik1)
  elpd_loo1 = sum(loos1)
  se_elpd_loo1 = std(loos1) * sqrt(n)
  @printf(">> elpd_loo = %.1f, SE(elpd_loo) = %.1f\n", elpd_loo1, se_elpd_loo1)

  # Check the shape parameter k of the generalized Pareto distribution
  if all(pk1 .< 0.5)
      println("All Pareto k estimates OK (k < 0.5).")
  else
    pk_good = sum(pk1 .<= 0.5)
    pk_ok = length(pk1[pk1 .<= 0.7]) - pk_good
    pk_bad = length(pk1[pk1 .<= 1]) - pk_good - pk_ok
    println((good=pk_good, ok=pk_ok, bad=pk_bad, very_bad=sum(pk1 .> 1)))
  end
end

begin
  scatter(pk1, xlab="Datapoint", ylab="Pareto shape k",
    marker=2.5, lab="Pk points")
  hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
  hline!([1], lab="pk = 1.0")
  title!("PSIS diagnostic plot for poisson-log model.")
  savefig(joinpath(ProjDir, "diag_plot_1.png"))
end

#=
*Simple negative binomial regression example
*using the 2nd parametrization of the negative
*binomial distribution, see section 40.1-3 in the Stan
*reference guide
=#

roaches2_stan = "
data {
  int<lower=0> N;
  int<lower=0> K;
  vector[N] exposure2;
  vector[N] roach1;
  vector[N] senior;
  vector[N] treatment;
  int y[N];
}
transformed data {
  vector[N] log_expo;
  log_expo = log(exposure2);
}
parameters {
  real phi;
  vector[K] beta;
}
transformed parameters {
   vector[N] eta;
   eta = log_expo + beta[1] + beta[2] * roach1 + beta[3] * treatment
                + beta[4] * senior;
}
model {
  phi ~ normal(0, 10);
  beta[1] ~ cauchy(0,10);   //prior for the intercept following Gelman 2008
  for(i in 2:K)
   beta[i] ~ cauchy(0,2.5); //prior for the slopes following Gelman 2008
  y ~ neg_binomial_2_log(eta, phi);
}
generated quantities {
 vector[N] log_lik;
 for(i in 1:N){
  log_lik[i] <- neg_binomial_2_log_lpmf(y[i] | eta[i], phi);
 }
}
";

tmpdir=joinpath(ProjDir, "tmp")
sm2 = SampleModel("roaches2", roaches2_stan; tmpdir)
rc2 = stan_sample(sm2; data)

if success(rc2)
  read_summary(sm2, true)
  nt2 = read_samples(sm2)

  # Compute LOO and standard error
  log_lik2 = nt2.log_lik'
  loo2, loos2, pk2 = psisloo(log_lik2)
  elpd_loo2 = sum(loos2)
  se_elpd_loo2 = std(loos2) * sqrt(n)
  @printf(">> elpd_loo = %.1f, SE(elpd_loo) = %.1f\n", elpd_loo2, se_elpd_loo2)

  # Check the shape parameter k of the generalized Pareto distribution
  pk_qualify(pk2) |> display

end

begin
  scatter(pk2, xlab="Datapoint", ylab="Pareto shape k",
    marker=2.5, lab="Pk points")
  hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
  hline!([1], lab="pk = 1.0")
  title!("PSIS diagnostic plot for neg-binomial model.")
  savefig(joinpath(ProjDir, "diag_plot_2.png"))
end


text:
__precompile__()
module GtkSourceWidget

export GtkSourceLanguage, GtkSourceLanguageManager, GtkSourceBuffer,
       GtkSourceView, GtkSourceCompletionItem,
       GtkSourceStyleSchemeManager, GtkSourceStyleScheme, GtkSourceMap,
       GtkSourceSearchContext, GtkSourceSearchSettings

export scheme, language, show_line_numbers!, auto_indent!, style_scheme, style_scheme!,
       max_undo_levels, max_undo_levels!, undo!, redo!, canundo, canredo, undomanager,
       highlight_current_line!, highlight_matching_brackets, source_view_get_gutter,
       reset_undomanager, set_view, get_view, style_scheme_chooser, style_scheme_chooser_button,
       set_search_text, get_search_text, search_context_forward, highlight, search_context_replace

using Gtk
using Compat

import ..Gtk: suffix, GObject, GtkTextIter

const MutableGtkTextIter = Gtk.GLib.MutableTypes.Mutable{GtkTextIter}

mutable(it::GtkTextIter) = Gtk.GLib.MutableTypes.mutable(it)

if Gtk.gtk_version == 3
    @static if is_windows()
        const libgtksourceview = Pkg.dir() * "/GtkSourceWidget/bin/libgtksourceview-3.0-1.dll"
    end
    @static if is_linux()
        try
            strip(readstring(pipeline(`ldconfig -p`, `grep libgtksourceview-3`, `cut -d'>' -f2`)))
        catch
            run(`sudo apt-get install libgtksourceview-3.0-1`)
        end
        const libgtksourceview = strip(readstring(pipeline(`ldconfig -p`, `grep libgtksourceview-3`, `cut -d'>' -f2`, `head -1`)))
	end
    @static if is_apple()
        if !isfile( Pkg.dir() * "/Homebrew/deps/usr/lib/libgtksourceview-3.0.dylib" )
            using Homebrew
            Homebrew.add("gtksourceview3")
        end
        const libgtksourceview = Pkg.dir() * "/Homebrew/deps/usr/lib/libgtksourceview-3.0"
    end
else
    error("Unsupported Gtk version $gtk_version")
end

### GtkSourceLanguage

Gtk.@Gtype GtkSourceLanguage libgtksourceview gtk_source_language

### GtkSourceLanguageManager

Gtk.@Gtype GtkSourceLanguageManager libgtksourceview gtk_source_language_manager
function GtkSourceLanguageManagerLeaf(default=true)
  if default
    GtkSourceLanguageManagerLeaf(
      ccall((:gtk_source_language_manager_get_default,libgtksourceview),Ptr{GObject},()))
  else
    GtkSourceLanguageManagerLeaf(
      ccall((:gtk_source_language_manager_new,libgtksourceview),Ptr{GObject},()))
  end
end

get_search_path(manager::GtkSourceLanguageManager) = ccall((:gtk_source_language_manager_get_search_path,libgtksourceview),Ptr{Ptr{UInt8}},
    (Ptr{GObject},),manager)

    # 3-element Array{Ptr{UInt8},1}:
    # Ptr{UInt8} @0x000000001bb85ac0
    # Ptr{UInt8} @0x000000001bb84e70
    # Ptr{UInt8} @0x0000000000000000
    #
    # julia> bytestring(pointer_to_array(p,2)[1])

set_search_path(manager::GtkSourceLanguageManager,dir)  = ccall((:gtk_source_language_manager_set_search_path,libgtksourceview),Void,
    (Ptr{GObject}, Ptr{Ptr{UInt8}}),manager, dir)

# Use: GtkSourceWidget.set_search_path(m,Any["c:/ad","yp",C_NULL])

language(manager::GtkSourceLanguageManager, id::String) = GtkSourceLanguage(
    ccall((:gtk_source_language_manager_get_language,libgtksourceview),Ptr{GObject},
          (Ptr{GObject},Ptr{UInt8}),manager,string(id)))

### GtkSourceStyle

Gtk.@Gtype GtkSourceStyle libgtksourceview gtk_source_style

### GtkSourceStyleScheme

Gtk.@Gtype GtkSourceStyleScheme libgtksourceview gtk_source_style_scheme

style(scheme::GtkSourceStyleScheme, style_id::String)  = GtkSourceStyle(
    ccall((:gtk_source_style_scheme_get_style,libgtksourceview),Ptr{GObject},
    (Ptr{GObject}, Ptr{UInt8}), scheme, string(style_id)))

### GtkSourceStyleSchemeManager

Gtk.@Gtype GtkSourceStyleSchemeManager libgtksourceview gtk_source_style_scheme_manager
function GtkSourceStyleSchemeManagerLeaf(default=true)
  if default
    GtkSourceStyleSchemeManagerLeaf(
      ccall((:gtk_source_style_scheme_manager_get_default,libgtksourceview),Ptr{GObject},()))
  else
    GtkSourceStyleSchemeManagerLeaf(
      ccall((:gtk_source_style_scheme_manager_get_new,libgtksourceview),Ptr{GObject},()))
  end
end

style_scheme(manager::GtkSourceStyleSchemeManager, scheme_id::String) = GtkSourceStyleScheme(
    ccall((:gtk_source_style_scheme_manager_get_scheme,libgtksourceview),Ptr{GObject},
          (Ptr{GObject},Ptr{UInt8}),manager,string(scheme_id)))

set_search_path(manager::GtkSourceStyleSchemeManager,dir)  = ccall((:gtk_source_style_scheme_manager_set_search_path,libgtksourceview),Void,
    (Ptr{GObject}, Ptr{Ptr{UInt8}}),manager, dir)

### GtkSourceUndoManager

### This is a hack!!!
type GtkSourceUndoManagerI
  handle::Ptr{Void}
end

undo!(manager::GtkSourceUndoManagerI) =
  ccall((:gtk_source_undo_manager_undo,libgtksourceview),Void,
        (Ptr{Void},),manager.handle)

redo!(manager::GtkSourceUndoManagerI) =
  ccall((:gtk_source_undo_manager_redo,libgtksourceview),Void,
        (Ptr{Void},),manager.handle)

canundo(manager::GtkSourceUndoManagerI) = Bool(
  ccall((:gtk_source_undo_manager_can_undo,libgtksourceview),Cint,
        (Ptr{Void},),manager.handle))

canredo(manager::GtkSourceUndoManagerI) = Bool(
  ccall((:gtk_source_undo_manager_can_redo,libgtksourceview),Cint,
        (Ptr{Void},),manager.handle))

### GtkSourceBuffer

Gtk.@Gtype GtkSourceBuffer libgtksourceview gtk_source_buffer
GtkSourceBufferLeaf() = GtkSourceBufferLeaf(
     ccall((:gtk_source_buffer_new,libgtksourceview),Ptr{GObject}, (Ptr{Void},), C_NULL) )

GtkSourceBufferLeaf(lang::GtkSourceLanguage) = GtkSourceBufferLeaf(
     ccall((:gtk_source_buffer_new_with_language,libgtksourceview),Ptr{GObject},
           (Ptr{GObject},), lang) )

highlight_syntax(buffer::GtkSourceBuffer, highlight::Bool) =
  ccall((:gtk_source_buffer_set_highlight_syntax,libgtksourceview),Ptr{GObject},
        (Ptr{GObject},Cint),buffer,Int32(highlight))

highlight_matching_brackets(buffer::GtkSourceBuffer, highlight::Bool) =
  ccall((:gtk_source_buffer_set_highlight_matching_brackets,libgtksourceview),Ptr{GObject},
        (Ptr{GObject},Cint),buffer,Int32(highlight))

style_scheme!(buffer::GtkSourceBuffer, scheme::GtkSourceStyleScheme) =
  ccall((:gtk_source_buffer_set_style_scheme,libgtksourceview),Void,
        (Ptr{GObject},Ptr{GObject}),buffer,scheme)

max_undo_levels(buffer::GtkSourceBuffer) =
  ccall((:gtk_source_buffer_get_max_undo_levels,libgtksourceview),Cint,
        (Ptr{GObject},),buffer)

max_undo_levels!(buffer::GtkSourceBuffer, levels::Integer) =
  ccall((:gtk_source_buffer_set_max_undo_levels,libgtksourceview),Void,
        (Ptr{GObject},Cint),buffer,levels)

undo!(buffer::GtkSourceBuffer) =
  ccall((:gtk_source_buffer_undo,libgtksourceview),Void,
        (Ptr{GObject},),buffer)

redo!(buffer::GtkSourceBuffer) =
  ccall((:gtk_source_buffer_redo,libgtksourceview),Void,
        (Ptr{GObject},),buffer)

canundo(buffer::GtkSourceBuffer) = Bool(
  ccall((:gtk_source_buffer_can_undo,libgtksourceview),Cint,
        (Ptr{GObject},),buffer))

canredo(buffer::GtkSourceBuffer) = Bool(
  ccall((:gtk_source_buffer_can_redo,libgtksourceview),Cint,
        (Ptr{GObject},),buffer))

undomanager(buffer::GtkSourceBuffer) = GtkSourceUndoManagerI(
  ccall((:gtk_source_buffer_get_undo_manager,libgtksourceview),Ptr{Void},
        (Ptr{GObject},),buffer))

reset_undomanager(buffer::GtkSourceBuffer) = ccall((:gtk_source_buffer_set_undo_manager,libgtksourceview),Void,
        (Ptr{GObject},Ptr{GObject}),buffer,C_NULL)

### GtkSourceView

Gtk.@Gtype GtkSourceView libgtksourceview gtk_source_view
GtkSourceViewLeaf(buffer=GtkSourceBuffer()) = GtkSourceView(
    ccall((:gtk_source_view_new_with_buffer,libgtksourceview),Ptr{GObject},(Ptr{GObject},),buffer))

show_line_numbers!(view::GtkSourceView, show::Bool) =
  ccall((:gtk_source_view_set_show_line_numbers,libgtksourceview),Void,
        (Ptr{GObject},Cint),view,Int32(show))

auto_indent!(view::GtkSourceView, auto::Bool) =
  ccall((:gtk_source_view_set_auto_indent,libgtksourceview),Void,
        (Ptr{GObject},Cint),view,Int32(auto))

highlight_current_line!(view::GtkSourceView, hl::Bool) =
  ccall((:gtk_source_view_set_highlight_current_line,libgtksourceview),Void,
        (Ptr{GObject},Cint),view,Int32(hl))

### GtkSourceCompletion

Gtk.@Gtype GtkSourceCompletion libgtksourceview gtk_source_completion

get_completion(view::GtkSourceView) =
    ccall((:gtk_source_view_get_completion,libgtksourceview),Ptr{GtkSourceCompletion},
        (Ptr{GObject},),view)

### GtkSourceCompletionContext

baremodule GtkSourceCompletionActivation
    const GTK_SOURCE_COMPLETION_ACTIVATION_NONE = 0
    const GTK_SOURCE_COMPLETION_ACTIVATION_INTERACTIVE = 1
    const GTK_SOURCE_COMPLETION_ACTIVATION_USER_REQUESTED = 2
end

Gtk.@Gtype GtkSourceCompletionContext libgtksourceview gtk_source_completion_context

#function add_proposals(context::GtkSourceCompletionContext, provider::GtkSourceCompletionProvider,
#                       proposals::GList, finished::Bool)
#  ccall((:gtk_source_completion_context_add_proposals,libgtksourceview),Void,
#        (Ptr{GObject},Ptr{GObject}, Ptr{GList},Cint),context, provider, &proposals, finished)
#end

#getiter(context::GtkSourceCompletionContext, iter::GtkTextIter) =
#   ccall((:gtk_source_completion_context_get_iter,libgtksourceview),Void,
#        (Ptr{GObject},Ptr{Void}),context, iter)

getactivation(context::GtkSourceCompletionContext) =
   ccall((:gtk_source_completion_context_get_activation,libgtksourceview),Cint,
        (Ptr{GObject},),context)

### GtkSourceCompletionItem

Gtk.@Gtype GtkSourceCompletionItem libgtksourceview gtk_source_completion_item

GtkSourceCompletionItemLeaf(label::String, text::String) = GtkSourceCompletionItemLeaf(
    ccall((:gtk_source_completion_item_new,libgtksourceview),Ptr{GObject},
        (Ptr{UInt8},Ptr{UInt8},Ptr{Void},Ptr{UInt8}),
        string(label),string(text),C_NULL,C_NULL))

### GtkSourceCompletionProvider

Gtk.GLib.@Giface GtkSourceCompletionProvider libgtksourceview gtk_source_completion_provider

## GtkSourceSearchSettings

Gtk.@Gtype GtkSourceSearchSettings libgtksourceview gtk_source_search_settings

GtkSourceSearchSettingsLeaf() = GtkSourceSearchSettingsLeaf(
    ccall((:gtk_source_search_settings_new,libgtksourceview),Ptr{GObject},())
)

set_search_text(settings::GtkSourceSearchSettings, text::String) =
    ccall((:gtk_source_search_settings_set_search_text,libgtksourceview),Void,
        (Ptr{GObject},Ptr{UInt8}),settings,text)

function get_search_text(settings::GtkSourceSearchSettings)

    s = ccall((:gtk_source_search_settings_get_search_text,libgtksourceview),Ptr{UInt8},
    (Ptr{GObject},),settings)

    return s == C_NULL ? "" : string(s)
end

## GtkSourceSearchContext

Gtk.@Gtype GtkSourceSearchContext libgtksourceview gtk_source_search_context

GtkSourceSearchContextLeaf(buffer::GtkSourceBuffer,settings::GtkSourceSearchSettings) = GtkSourceSearchContextLeaf(
    ccall((:gtk_source_search_context_new,libgtksourceview),Ptr{GObject},(Ptr{GObject},Ptr{GObject}),buffer,settings))

search_context_forward(search::GtkSourceSearchContext, iter::GtkTextIter,match_start::MutableGtkTextIter,match_end::MutableGtkTextIter) =
    convert(Bool,ccall((:gtk_source_search_context_forward,libgtksourceview),Cint,
        (Ptr{GObject},Ref{GtkTextIter},Ptr{MutableGtkTextIter},Ptr{MutableGtkTextIter}),
        search,iter,match_start,match_end))

function search_context_forward(search::GtkSourceSearchContext, iter::GtkTextIter)

    match_start = mutable(GtkTextIter())
    match_end   = mutable(GtkTextIter())

    found = search_context_forward(search,iter,match_start,match_end)
    return (found,match_start,match_end)
end

function search_context_replace(
    search::GtkSourceSearchContext,
    match_start::MutableGtkTextIter, match_end::MutableGtkTextIter,
    replace::String)

    out = ccall((:gtk_source_search_context_replace,libgtksourceview),Cint,
        (Ptr{GObject},Ptr{MutableGtkTextIter},Ptr{MutableGtkTextIter},Ptr{UInt8},Cint,Ptr{Void}),
        search,match_start,match_end,string(replace),-1,C_NULL)

    return convert(Bool,out)
end

function search_context_replace_all(search::GtkSourceSearchContext, replace::String)

    out = ccall((:gtk_source_search_context_replace_all,libgtksourceview),Cint,
        (Ptr{GObject},Ptr{UInt8},Cint,Ptr{Void}),
        search,string(replace),-1,C_NULL)

    return out
end

highlight(search::GtkSourceSearchContext, highlight::Bool) =
    ccall((:gtk_source_search_context_set_highlight,libgtksourceview),Void,
        (Ptr{GObject},Cint), search,highlight)

### GtkSourceGutter

source_view_get_gutter(view::GtkSourceView) =  ccall((:gtk_source_view_get_gutter,libgtksourceview),Ptr{Void},(Ptr{GObject},Cint),view,Int32(1))

### GtkSourceMap

function _define_source_map()
    if is_linux()
        return :(
            global const SOURCE_MAP = false
        )
    else
        return quote
            global const SOURCE_MAP  = true
            Gtk.@Gtype GtkSourceMap libgtksourceview gtk_source_map

    		GtkSourceMapLeaf() = GtkSourceMapLeaf( ccall((:gtk_source_map_new,libgtksourceview),Ptr{GObject},()) )

    		get_view(map::GtkSourceMap) =
    		    ccall((:gtk_source_map_get_view,libgtksourceview),Ptr{GtkSourceView},
    		        (Ptr{GObject},),map)

    		set_view(map::GtkSourceMap,view::GtkSourceView) =
    		    ccall((:gtk_source_map_set_view,libgtksourceview),Void,
    		        (Ptr{GObject},Ptr{GObject}),map,view)
        end
    end
end

macro define_source_map() esc(_define_source_map()) end

@define_source_map()

## GtkSourceStyleSchemeChooserWidget

style_scheme_chooser() = Gtk.GtkWidgetLeaf( ccall((:gtk_source_style_scheme_chooser_widget_new,libgtksourceview),Ptr{GObject},()) )
style_scheme_chooser_button() = Gtk.GtkWidgetLeaf( ccall((:gtk_source_style_scheme_chooser_button_new,libgtksourceview),Ptr{GObject},()) )

style_scheme(chooser::Gtk.GtkWidget) = GtkSourceStyleScheme(
    ccall((:gtk_source_style_scheme_chooser_get_style_scheme,libgtksourceview),Ptr{GObject},
        (Ptr{GObject},),chooser)  )

#? https://github.com/JuliaLang/Gtk.jl/blob/40f7442dea20919c5a7b137594fb6d7636fa2329/src/selectors.jl

function __init__()

    global const sourceLanguageManager = GtkSourceLanguageManager()
    GtkSourceWidget.set_search_path(sourceLanguageManager,
      Any[Pkg.dir() * "/GtkSourceWidget/share/gtksourceview-3.0/language-specs/",C_NULL])
end

end # module


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

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

# ╔═╡ ff275f7e-b9fe-11eb-1f20-39ca812b0fbe
using WAV, PlutoUI, DSP, Plots,UnicodePlots

# ╔═╡ 983dde83-0ad7-4d1f-81b8-a85e74b137f0
gr()

# ╔═╡ 59209394-0315-4e4e-bd83-8b8d08b9acd4
begin
# # fs = 8e3
# t = 0.0:1/fs:prevfloat(1.0)
# f = 1e3
# y = sin.(2pi * f * t) * 0.1
# wavwrite(y, "example.wav", Fs=fs)

# y, fs = wavread("example.wav")
# y = sin.(2pi * 2f * t) * 0.1
# wavappend(y, "example.wav")

y, fs1 = wavread("piano-phrase.wav")
wavplay(y, fs1)
end

# ╔═╡ 8413a04a-a377-40e1-b798-26d8964f8ada
y #we use 64-bit float

# ╔═╡ 9ea99bc6-38d0-40bf-b498-5aff2cf957d1
fs1 #sampling rate

# ╔═╡ 13937baf-4bb7-4de8-9a8d-9a6d0228fb52
size(y)

# ╔═╡ da16ccb9-b5a8-4875-94d8-9c60538ae41c

begin
plot(0:1/fs1:(length(y)-1)/fs1, y)
# xlabel("Time [s]")
end

# ╔═╡ 4c7ab8a6-8fbb-4024-b434-ab512d156f61
wavplay(y, fs1)

# ╔═╡ fecbcf72-4146-4dff-a2f7-6e5c5e07bfe1
@bind A Slider(.1:.1:1,default=.8,show_value=true)

# ╔═╡ 6a68a83a-6f4c-4355-8ddc-797cccaac1f4
@bind phi Slider(0:pi/6:2pi,default=pi/2,show_value=true)

# ╔═╡ a3cbf87e-ae75-4418-b5b4-44d6ac5f1469
# @bind f0 Slider(0:10:1000,default=440,show_value=true)

# ╔═╡ f4b75e73-47f9-4786-8210-8619daf50cb5
@bind fsd Slider(0:1000:50000,default=44100,show_value=true)

# ╔═╡ e7426fcc-bb42-4df0-8390-ea2085021860
begin
f0 = 880
t = 1:1/fsd:11
wave = A./t .*sin.(2*pi*f0*t .+ phi)

#y1 = squarewave.(f1*t1) * 1
#y1 = (sin.(2pi * f1* t1) .*0.1)
#y1 = trianglewave.(f1* t1) * 1
# y1 = (sin.(2pi * f1* t1) .*t1) .+ (sin.(2pi * f1*2* t1 .+pi) .*t1/2)
wavplay(wave, fsd)
end

# ╔═╡ 521e87e1-d9eb-4bfc-96f5-b8678c7c83a1
lineplot(t,wave)

# ╔═╡ 9b852298-1cb9-49d9-9384-1be059d2918a
begin
@bind f1 Slider(0:10:1000,default=200,show_value=true)
end

# ╔═╡ 24ecd45a-9863-4cc2-b23c-ceb2c2f8c9e9
@bind f2 Slider(0:10:1000,default=400,show_value=true)

# ╔═╡ 5398de30-c98f-405f-a34c-a790a6458e17
@bind f3 Slider(0:10:1000,default=600,show_value=true)

# ╔═╡ 639ac6cb-4763-406e-8f31-8c190ef2ea30
@bind f4 Slider(0:10:1000,default=800,show_value=true)

# ╔═╡ 77737e79-de60-4134-a95e-2e40b2ca1ec3
@bind f5 Slider(0:10:1000,default=1000,show_value=true)

# ╔═╡ 58c40dd6-23f6-4a97-8389-acff4c906a4c
begin
ta = 0:1/fsd:0.008
harmonic1 = A*sin.(2*pi*f1*ta .+ phi)
harmonic2 = A*sin.(2*pi*f2*ta .+ phi)
harmonic3 = A*sin.(2*pi*f3*ta .+ phi)
harmonic4 = A*sin.(2*pi*f4*ta .+ phi)
harmonic5 = A*sin.(2*pi*f5*ta .+ phi)
end

# ╔═╡ 537e7423-cb27-49b1-8cec-947051ba2964
lineplot(ta,harmonic1)

# ╔═╡ c25eb43b-2fc7-4f96-a4f4-f48f89138786
lineplot(ta,harmonic2)

# ╔═╡ 4bf9ddc2-c3d9-48af-9afd-509aa0505e53
lineplot(ta,harmonic3)

# ╔═╡ cabb81ec-1d67-46f4-b4e3-296f763fc644
lineplot(ta,harmonic4)

# ╔═╡ 0bc6ab88-9498-417e-95af-4b733a0db37a
lineplot(ta,harmonic5)

# ╔═╡ 938ed7b9-2a5e-4a8a-a8fe-f3e1e179324e
composite = harmonic1 .+ harmonic2 .+ harmonic3 .+ harmonic4 .+harmonic5

# ╔═╡ cb7abc8f-259d-4509-8e67-69dc7fbdc2bc
lineplot(ta,composite)

# ╔═╡ Cell order:
# ╠═ff275f7e-b9fe-11eb-1f20-39ca812b0fbe
# ╠═983dde83-0ad7-4d1f-81b8-a85e74b137f0
# ╠═59209394-0315-4e4e-bd83-8b8d08b9acd4
# ╠═8413a04a-a377-40e1-b798-26d8964f8ada
# ╠═9ea99bc6-38d0-40bf-b498-5aff2cf957d1
# ╠═13937baf-4bb7-4de8-9a8d-9a6d0228fb52
# ╠═da16ccb9-b5a8-4875-94d8-9c60538ae41c
# ╠═4c7ab8a6-8fbb-4024-b434-ab512d156f61
# ╠═fecbcf72-4146-4dff-a2f7-6e5c5e07bfe1
# ╠═6a68a83a-6f4c-4355-8ddc-797cccaac1f4
# ╠═a3cbf87e-ae75-4418-b5b4-44d6ac5f1469
# ╠═f4b75e73-47f9-4786-8210-8619daf50cb5
# ╠═e7426fcc-bb42-4df0-8390-ea2085021860
# ╠═521e87e1-d9eb-4bfc-96f5-b8678c7c83a1
# ╠═9b852298-1cb9-49d9-9384-1be059d2918a
# ╠═24ecd45a-9863-4cc2-b23c-ceb2c2f8c9e9
# ╠═5398de30-c98f-405f-a34c-a790a6458e17
# ╠═639ac6cb-4763-406e-8f31-8c190ef2ea30
# ╠═77737e79-de60-4134-a95e-2e40b2ca1ec3
# ╠═58c40dd6-23f6-4a97-8389-acff4c906a4c
# ╠═537e7423-cb27-49b1-8cec-947051ba2964
# ╠═c25eb43b-2fc7-4f96-a4f4-f48f89138786
# ╠═4bf9ddc2-c3d9-48af-9afd-509aa0505e53
# ╠═cabb81ec-1d67-46f4-b4e3-296f763fc644
# ╠═0bc6ab88-9498-417e-95af-4b733a0db37a
# ╠═938ed7b9-2a5e-4a8a-a8fe-f3e1e179324e
# ╠═cb7abc8f-259d-4509-8e67-69dc7fbdc2bc

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
module CourseMatch

	VERSION = VersionNumber(0,0,2)

using JuMP, Cbc, CSV, DataFrames, SparseArrays, Dates
using DataStructures, Distributions, DataFramesMeta
using JuMP, Cbc

# includes
include("student.jl")
include("demand.jl")

include("clearing.jl")

# exports

end # module


text:
module JudiLing

using DataFrames
using Random, Distributions
using SparseArrays, LinearAlgebra, Statistics, SuiteSparse
using BenchmarkTools
using DataStructures
using ProgressMeter
using CSV
using GZip
using PyCall

include("utils.jl")
include("pyndl.jl")
include("wh.jl")
include("make_cue_matrix.jl")
include("make_semantic_matrix.jl")
include("cholesky.jl")
include("make_adjacency_matrix.jl")
include("make_yt_matrix.jl")
include("find_path.jl")
include("eval.jl")
include("output.jl")
include("preprocess.jl")
include("pickle.jl")
include("test_combo.jl")
include("display.jl")

end