text:
using Plots
using DataStructures
using Statistics
using LaTeXStrings
L=200
function BoundaryCondition(i)
    if i==L+1
        return 1
    elseif i==0
        return L
    else
        return i
    end
end
function height_cal(i,arr) #calculating the height of the layer in a column
    row=1
    height=0
    for j in 1:L
        if arr[row,i]!=0.0
            height+=1
            row+=1
        else
            break #exit the loop for the next column
        end
    end
    return height
end
function ColumnChoosing(arr, col) #choosing which column to fall into
    h_dict=DataStructures.OrderedDict() #an ordered dictionary : height=>column_num
    for i in (col-1):1:(col+1)
        i=BoundaryCondition(i)
        height= height_cal(i,arr)
        h_dict[height]=i
    end
    hArr=collect(keys(h_dict)) #having all the heights in an array
    return get(h_dict,min(hArr...), "ERROR") #returning the column respectively
end
function std_calculater(arr,stdList)
    hight=[]
    for col in eachcol(arr)
        count=0
        for element in col
            if element!=0
                count+=1
            end
        end
        push!(hight, count)
    end
    std_num=std(hight)
    push!(stdList, std_num)
    return stdList
end
L=200
arr=zeros((L,L))
color=1
stdList=[]
t_interval=ceil.(Int, exp.(1:0.5:13))
for i in 1:24
    falling_n=t_interval[i+1]-t_interval[i]
    for each in 1:falling_n
        column=rand(1:L)
        col=ColumnChoosing(arr, column)
        row=1
        for j in 1:L
            if arr[row,col]==0.0
                arr[row,col]=color
                break
            else
                row+=1
            end
        end
        if i%(10*200*color)==0
            color+=1
        end
    end
    stdList=std_calculater(arr,stdList)
end
BigStdList=[]
run_number=100
for num in 1:run_number
    println("run number=", num)
    L=200
    arr=zeros((L,L))
    color=1
    stdList=[]
    time=ceil.(Int, exp.(0:0.5:13))
    for i in 1:24
        falling_n=time[i+1]-time[i]
        for each in 1:falling_n
            column=rand(1:L)
            col=ColumnChoosing(arr, column)
            row=1
            for j in 1:L
                if arr[row,col]==0.0
                    arr[row,col]=color
                    break
                else
                    row+=1
                end
            end
            if i%(10*200*color)==0
                color+=1
            end
        end
        stdList=std_calculater(arr,stdList)
    end
    push!(BigStdList, stdList)
end
BiggerMeanList=[]
BiggerStdList=[]
for i in 1:24
    tempvalue=[]
    for j in 1:100
        push!(tempvalue,BigStdList[j][i])
    end
    push!(BiggerStdList, std(tempvalue))
    push!(BiggerMeanList, mean(tempvalue))
end
function LineFit(Time, Mean)
    a=[hcat(log.(Time)) reshape(ones(13), 13, 1)]
    b=hcat(log.(Mean))
    line=(a\b)
    x = 1:10
    y = x .* line[1] .+ line[2]
    return x, y, line
end
BiggerMeanList[7:19]
X,Y,Line=LineFit(t_interval[7:19], BiggerMeanList[7:19])
plot(log.(t_interval)[1:24], log.(BiggerMeanList)[1:24],yerr=BiggerStdList, xlabel="time", ylabel="w(t)")
plot!(X,Y,c= :black,label = L"y = %$(round(Line[1],digits= 2))x + %$(round(Line[2],digits= 2))", legend=:bottomright)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\RandomBallisticDepositionWithRelaxation3.png")


text:
###
#   IntervalUnion Arithmetic
###

# Special case for /. Probably better with promotion rules and convert
for op in (:+, :-, :/, :*, :min, :max, :^, :log, :<, :>)

    if op != :/; @eval ($op)( x::IntervalU, y::IntervalU) = intervalU([$op(xv, yv) for xv in x.v, yv in y.v][:]); end

    @eval ($op)( x::IntervalU, y::Interval) = intervalU(broadcast($op, x.v, y))
    if op != :/; @eval ($op)( x::Interval, y::IntervalU) = intervalU(broadcast($op, x, y.v)); end

    @eval ($op)( x::IntervalU, n::Real) = intervalU(broadcast($op, x.v, n))
    @eval ($op)( n::Real, x::IntervalU) = intervalU(broadcast($op, n, x.v))

end

for op in (:-, :sin, :cos, :tan, :exp, :log)
    @eval ($op)( x::IntervalU) = intervalU(broadcast($op, x.v))
end

# Does x/y for y\{0} if 0 ∈ y
function /(x :: IntervalU, y :: IntervalU)
    yNew = bisect(y, 0)
    return intervalU(broadcast(/, x.v, yNew.v))
end

function /(x :: Interval, y :: IntervalU)
    yNew = bisect(y, 0)
    return intervalU(broadcast(/, x, yNew.v))
end

function sqrt( x:: IntervalU)
    us = sqrt.(x.v)
    return intervalU([us; -us])
end

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to IAU-2006 transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/iau2006/iau2006-cio.jl
# ==================================================

# Functions r_itrf_to_tirs_iau2006 and r_tirs_to_itrf_iau2006
# -----------------------------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   JD_TT  = 2453101.828154745
#   x_p    = -0.140682"
#   y_p    = +0.333309"
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
# one gets the following:
#
#   r_tirs  = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_tirs  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function r_itrf_to_tirs_iau2006 and r_tirs_to_itrf_iau2006" begin
    JD_TT = 2453101.828154745
    x_p   = -0.140682*pi/(180*3600)
    y_p   = +0.333309*pi/(180*3600)

    ## r_itrf_to_tirs_iau2006
    ## ===================

    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]
    v_itrf = [-3.225636520; -2.872451450; +5.531924446]

    ## DCM
    ## ---

    D_TIRS_ITRF = r_itrf_to_tirs_iau2006(JD_TT, x_p, y_p)

    r_tirs = D_TIRS_ITRF*r_itrf
    v_tirs = D_TIRS_ITRF*v_itrf

    @test r_tirs[1] ≈ -1033.47503120 atol=1e-7
    @test r_tirs[2] ≈ +7901.30558560 atol=1e-7
    @test r_tirs[3] ≈ +6380.34453270 atol=1e-7

    @test v_tirs[1] ≈ -3.2256327470  atol=1e-9
    @test v_tirs[2] ≈ -2.8724425110  atol=1e-9
    @test v_tirs[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TIRS_ITRF = r_itrf_to_tirs_iau2006(Quaternion, JD_TT, x_p, y_p)

    r_tirs = vect(q_TIRS_ITRF \ r_itrf * q_TIRS_ITRF)
    v_tirs = vect(q_TIRS_ITRF \ v_itrf * q_TIRS_ITRF)

    @test r_tirs[1] ≈ -1033.47503120 atol=1e-7
    @test r_tirs[2] ≈ +7901.30558560 atol=1e-7
    @test r_tirs[3] ≈ +6380.34453270 atol=1e-7

    @test v_tirs[1] ≈ -3.2256327470  atol=1e-9
    @test v_tirs[2] ≈ -2.8724425110  atol=1e-9
    @test v_tirs[3] ≈ +5.5319312880  atol=1e-9

    ## r_tirs_to_itrf_iau2006
    ## ==============

    r_tirs  = [-1033.47503120; 7901.30558560; 6380.34453270]
    v_tirs  = [-3.2256327470; -2.8724425110; +5.5319312880]

    ## DCM
    ## ---

    D_ITRF_TIRS = r_tirs_to_itrf_iau2006(JD_TT, x_p, y_p)
    r_itrf = D_ITRF_TIRS*r_tirs
    v_itrf = D_ITRF_TIRS*v_tirs

    @test r_itrf[1] ≈ -1033.4793830 atol=1e-7
    @test r_itrf[2] ≈ +7901.2952754 atol=1e-7
    @test r_itrf[3] ≈ +6380.3565958 atol=1e-7

    @test v_itrf[1] ≈ -3.225636520  atol=1e-9
    @test v_itrf[2] ≈ -2.872451450  atol=1e-9
    @test v_itrf[3] ≈ +5.531924446  atol=1e-9

    ## Quaternion
    ## ----------

    q_ITRF_TIRS = r_tirs_to_itrf_iau2006(Quaternion, JD_TT, x_p, y_p)
    r_itrf = vect(q_ITRF_TIRS\r_tirs*q_ITRF_TIRS)
    v_itrf = vect(q_ITRF_TIRS\v_tirs*q_ITRF_TIRS)

    @test r_itrf[1] ≈ -1033.4793830 atol=1e-7
    @test r_itrf[2] ≈ +7901.2952754 atol=1e-7
    @test r_itrf[3] ≈ +6380.3565958 atol=1e-7

    @test v_itrf[1] ≈ -3.225636520  atol=1e-9
    @test v_itrf[2] ≈ -2.872451450  atol=1e-9
    @test v_itrf[3] ≈ +5.531924446  atol=1e-9
end

# Functions r_tirs_to_cirs_iau2006 and r_cirs_to_tirs_iau2006
# -----------------------------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   JD_UT1 = 2453101.827406783
#   LOD    = 0.0015563 s
#   r_tirs = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_tirs =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
# one gets the following:
#
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#   v_cirs =    -4.7453803300 i -    0.7903414530 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Functions r_tirs_to_cirs_iau2006 and r_cirs_to_tirs_iau2006" begin
    JD_UT1 = 2453101.827406783
    LOD    = 0.0015563
    w      = 7.292115146706979e-5*(1-LOD/86400)

    # r_tirs_to_cirs_iau2006
    # ===================

    r_tirs  = [-1033.47503120; 7901.30558560; 6380.34453270]
    v_tirs  = [-3.2256327470; -2.8724425110; +5.5319312880]

    # DCM
    # ---

    D_CIRS_TIRS = r_tirs_to_cirs_iau2006(JD_UT1)

    r_cirs = D_CIRS_TIRS * r_tirs
    v_cirs = D_CIRS_TIRS * (v_tirs + [0, 0, w] × r_tirs)

    @test r_cirs[1] ≈ +5100.01840470 atol=1e-7
    @test r_cirs[2] ≈ +6122.78636480 atol=1e-7
    @test r_cirs[3] ≈ +6380.34453270 atol=1e-7

    @test v_cirs[1] ≈ -4.7453803300  atol=1e-9
    @test v_cirs[2] ≈ +0.7903414530  atol=1e-9
    @test v_cirs[3] ≈ +5.5319312880  atol=1e-9

    # Quaternion
    # ----------

    q_CIRS_TIRS = r_tirs_to_cirs_iau2006(Quaternion, JD_UT1)

    r_cirs = vect(q_CIRS_TIRS\r_tirs * q_CIRS_TIRS)
    v_cirs = vect(q_CIRS_TIRS\(v_tirs + [0, 0, w] × r_tirs) * q_CIRS_TIRS)

    @test r_cirs[1] ≈ +5100.01840470 atol=1e-7
    @test r_cirs[2] ≈ +6122.78636480 atol=1e-7
    @test r_cirs[3] ≈ +6380.34453270 atol=1e-7

    @test v_cirs[1] ≈ -4.7453803300  atol=1e-9
    @test v_cirs[2] ≈ +0.7903414530  atol=1e-9
    @test v_cirs[3] ≈ +5.5319312880  atol=1e-9

    # r_cirs_to_tirs_iau2006
    # ===================

    r_cirs = [+5100.01840470; +6122.78636480; +6380.34453270]
    v_cirs = [-4.7453803300; +0.7903414530; +5.5319312880]

    # DCM
    # ---

    D_TIRS_CIRS = r_cirs_to_tirs_iau2006(JD_UT1)

    r_tirs = D_TIRS_CIRS * r_cirs
    v_tirs = D_TIRS_CIRS * v_cirs - [0, 0, w] × r_tirs

    @test r_tirs[1] ≈ -1033.47503120 atol=1e-7
    @test r_tirs[2] ≈ +7901.30558560 atol=1e-7
    @test r_tirs[3] ≈ +6380.34453270 atol=1e-7

    @test v_tirs[1] ≈ -3.2256327470  atol=1e-9
    @test v_tirs[2] ≈ -2.8724425110  atol=1e-9
    @test v_tirs[3] ≈ +5.5319312880  atol=1e-9

    # Quaternion
    # ----------

    q_TIRS_CIRS = r_cirs_to_tirs_iau2006(Quaternion, JD_UT1)

    r_tirs = vect(q_TIRS_CIRS \ r_cirs * q_TIRS_CIRS)
    v_tirs = vect(q_TIRS_CIRS \ v_cirs * q_TIRS_CIRS) - [0, 0, w] × r_tirs

    @test r_tirs[1] ≈ -1033.47503120 atol=1e-7
    @test r_tirs[2] ≈ +7901.30558560 atol=1e-7
    @test r_tirs[3] ≈ +6380.34453270 atol=1e-7

    @test v_tirs[1] ≈ -3.2256327470  atol=1e-9
    @test v_tirs[2] ≈ -2.8724425110  atol=1e-9
    @test v_tirs[3] ≈ +5.5319312880  atol=1e-9
end
#
# Functions r_cirs_to_gcrf_iau2006 and r_gcrf_to_cirs_iau2006
# -----------------------------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
# NOTE: It seems that the results in the Example 3-14 is computed **without**
# the `dX` and `dY` corrections, whereas in the Table 3-6 they are computed
# **with** the corrections.
#
#   JD_UT1 = 2453101.827406783
#   dX     = -0.000205"
#   dY     = -0.000136"
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#   v_cirs =    -4.7453803300 i -    0.7903414530 j +    5.5319312880 k [km/s]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#   v_gcrf =  -4.7432201610 i + 0.7905364950  j + 5.5337557240  k [km/s]
#
################################################################################

@testset "Functions r_cirs_to_gcrf_iau2006 and r_gcrf_to_cirs_iau2006" begin
    JD_TT = 2453101.827406783
    dX    = -0.000205 * pi / 180 / 3600
    dY    = -0.000136 * pi / 180 / 3600

    # r_cirs_to_gcrf_iau2006
    # ===================

    r_cirs = [+5100.01840470; +6122.78636480; +6380.34453270]
    v_cirs = [-4.7453803300; +0.7903414530; +5.5319312880]

    # DCM
    # ---

    D_GCRF_CIRS = r_cirs_to_gcrf_iau2006(JD_TT, dX, dY)
    r_gcrf = D_GCRF_CIRS * r_cirs
    v_gcrf = D_GCRF_CIRS * v_cirs

    @test r_gcrf[1] ≈ +5102.50895290 atol=9e-7
    @test r_gcrf[2] ≈ +6123.01139910 atol=9e-7
    @test r_gcrf[3] ≈ +6378.13693380 atol=9e-7

    @test v_gcrf[1] ≈ -4.7432201610  atol=1e-9
    @test v_gcrf[2] ≈ +0.7905364950  atol=1e-9
    @test v_gcrf[3] ≈ +5.5337557240  atol=1e-9

    # Quaternion
    # ----------

    q_GCRF_CIRS = r_cirs_to_gcrf_iau2006(Quaternion, JD_TT, dX, dY)

    r_gcrf = vect(q_GCRF_CIRS \ r_cirs * q_GCRF_CIRS)
    v_gcrf = vect(q_GCRF_CIRS \ v_cirs * q_GCRF_CIRS)

    @test r_gcrf[1] ≈ +5102.50895290 atol=9e-7
    @test r_gcrf[2] ≈ +6123.01139910 atol=9e-7
    @test r_gcrf[3] ≈ +6378.13693380 atol=9e-7

    @test v_gcrf[1] ≈ -4.7432201610  atol=1e-9
    @test v_gcrf[2] ≈ +0.7905364950  atol=1e-9
    @test v_gcrf[3] ≈ +5.5337557240  atol=1e-9

    # Function r_gcrf_to_cirs_iau2006
    # ============================

    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]
    v_gcrf = [-4.7432201610; 0.7905364950; 5.5337557240]

    # DCM
    # ---

    D_CIRS_GCRF = r_gcrf_to_cirs_iau2006(JD_TT, dX, dY)
    r_cirs = D_CIRS_GCRF * r_gcrf
    v_cirs = D_CIRS_GCRF * v_gcrf

    @test r_cirs[1] ≈ +5100.01840470 atol=9e-7
    @test r_cirs[2] ≈ +6122.78636480 atol=9e-7
    @test r_cirs[3] ≈ +6380.34453270 atol=9e-7

    @test v_cirs[1] ≈ -4.7453803300  atol=1e-9
    @test v_cirs[2] ≈ +0.7903414530  atol=1e-9
    @test v_cirs[3] ≈ +5.5319312880  atol=1e-9

    # Quaternion
    # ----------

    q_CIRS_GCRF = r_gcrf_to_cirs_iau2006(Quaternion, JD_TT, dX, dY)

    r_cirs = vect(q_CIRS_GCRF \ r_gcrf * q_CIRS_GCRF)
    v_cirs = vect(q_CIRS_GCRF \ v_gcrf * q_CIRS_GCRF)

    @test r_cirs[1] ≈ +5100.01840470 atol=9e-7
    @test r_cirs[2] ≈ +6122.78636480 atol=9e-7
    @test r_cirs[3] ≈ +6380.34453270 atol=9e-7

    @test v_cirs[1] ≈ -4.7453803300  atol=1e-9
    @test v_cirs[2] ≈ +0.7903414530  atol=1e-9
    @test v_cirs[3] ≈ +5.5319312880  atol=1e-9
end


text:
function SearchInterval(arg0::jdouble, arg1::jdouble)
    return SearchInterval((jdouble, jdouble), arg0, arg1)
end

function SearchInterval(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return SearchInterval((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_max(obj::SearchInterval)
    return jcall(obj, "getMax", jdouble, ())
end

function get_min(obj::SearchInterval)
    return jcall(obj, "getMin", jdouble, ())
end

function get_start_value(obj::SearchInterval)
    return jcall(obj, "getStartValue", jdouble, ())
end


text:
# This file was generated, do not modify it. # hide
xx = (Horsepower=range(50, 225, length=100) |> collect, )
yy = predict(mlm, xx)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy, lw=3)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g2.svg")) # hide


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SCIP_PaPILO_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SCIP_PaPILO")
JLLWrappers.@generate_main_file("SCIP_PaPILO", UUID("fc9abe76-a5e6-5fed-b0b7-a12f309cf031"))
end  # module SCIP_PaPILO_jll


text:
"""
    FeTerm{T,S}

Term with an explicit, constant matrix representation

Typically, an `FeTerm` represents the model matrix for the fixed effects.

!!! note
    `FeTerm` is not the same as [`FeMat`](@ref)!

# Fields
* `x`: full model matrix
* `piv`: pivot `Vector{Int}` for moving linearly dependent columns to the right
* `rank`: computational rank of `x`
* `cnames`: vector of column names
"""
struct FeTerm{T,S<:AbstractMatrix}
    x::S
    piv::Vector{Int}
    rank::Int
    cnames::Vector{String}
end

"""
    FeTerm(X::AbstractMatrix, cnms)

Convenience constructor for [`FeTerm`](@ref) that computes the rank and pivot with unit weights.

See the vignette "[Rank deficiency in mixed-effects models](@ref)" for more information on the
computation of the rank and pivot.
"""
function FeTerm(X::AbstractMatrix{T}, cnms) where {T}
    if iszero(size(X, 2))
        return FeTerm{T,typeof(X)}(X, Int[], 0, cnms)
    end
    rank, pivot = statsrank(X)
    # single-column rank deficiency is the result of a constant column vector
    # this generally happens when constructing a dummy response, so we don't
    # warn.
    if rank < length(pivot) && size(X, 2) > 1
        @warn "Fixed-effects matrix is rank deficient"
    end
    return FeTerm{T,typeof(X)}(X[:, pivot], pivot, rank, cnms[pivot])
end

"""
    FeTerm(X::SparseMatrixCSC, cnms)

Convenience constructor for a sparse [`FeTerm`](@ref) assuming full rank, identity pivot and unit weights.

Note: automatic rank deficiency handling may be added to this method in the future, as discused in
the vignette "[Rank deficiency in mixed-effects models](@ref)" for general `FeTerm`.
"""
function FeTerm(X::SparseMatrixCSC, cnms::AbstractVector{String})
    @debug "Full rank is assumed for sparse fixed-effect matrices."
    rank = size(X, 2)
    return FeTerm{eltype(X),typeof(X)}(X, collect(1:rank), rank, collect(cnms))
end

Base.copyto!(A::FeTerm{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.x, src)

Base.eltype(::FeTerm{T}) where {T} = T

function fullrankx(A::FeTerm)
    x, rnk = A.x, A.rank
    return rnk == size(x, 2) ? x : view(x, :, 1:rnk)  # this handles the zero-columns case
end

LinearAlgebra.rank(A::FeTerm) = A.rank

"""
    isfullrank(A::FeTerm)

Does `A` have full column rank?
"""
isfullrank(A::FeTerm) = A.rank == length(A.piv)

"""
    FeMat{T,S}

A matrix and a (possibly) weighted copy of itself.

Typically, an `FeMat` represents the fixed-effects model matrix with the response (`y`) concatenated as a final column.

!!! note
    `FeMat` is not the same as [`FeTerm`](@ref).

# Fields
- `xy`: original matrix, called `xy` b/c in practice this is `hcat(fullrank(X), y)`
- `wtxy`: (possibly) weighted copy of `xy` (shares storage with `xy` until weights are applied)

Upon construction the `xy` and `wtxy` fields refer to the same matrix
"""
mutable struct FeMat{T,S<:AbstractMatrix} <: AbstractMatrix{T}
    xy::S
    wtxy::S
end

function FeMat(A::FeTerm{T}, y::AbstractVector{T}) where {T}
    xy = hcat(fullrankx(A), y)
    return FeMat{T,typeof(xy)}(xy, xy)
end

Base.adjoint(A::FeMat) = Adjoint(A)

Base.eltype(::FeMat{T}) where {T} = T

Base.getindex(A::FeMat, i, j) = getindex(A.xy, i, j)

Base.length(A::FeMat) = length(A.xy)

function *(adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T}
    return adjoint(adjA.parent.wtxy) * B.wtxy
end

function LinearAlgebra.mul!(
    R::StridedVecOrMat{T}, A::FeMat{T}, B::StridedVecOrMat{T}
) where {T}
    return mul!(R, A.wtxy, B)
end

function LinearAlgebra.mul!(C, adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T}
    return mul!(C, adjoint(adjA.parent.wtxy), B.wtxy)
end

function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}
    if !isempty(sqrtwts)
        if A.xy === A.wtxy
            A.wtxy = similar(A.xy)
        end
        mul!(A.wtxy, Diagonal(sqrtwts), A.xy)
    end
    return A
end

Base.size(A::FeMat) = size(A.xy)

Base.size(A::FeMat, i::Integer) = size(A.xy, i)


text:
#
# Parrot.jl/docs/make.jl
#
# By Elijah Mathews (me@elijahmathews.com)
# openpgp4fpr:ac1d3fb1e8a5eb7d14bd587b2932c725055a90d8
#

push!(LOAD_PATH, "../src/")

using Documenter, Parrot

makedocs(
    format = Documenter.HTML(
        assets = ["assets/favicon.ico"],
    ),
    sitename = "Parrot",
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/elijahmathews/Parrot.jl.git",
    devbranch = "primary",
)


text:
# See "Statistical Applications of the multivariate skew-normal distribution" by Azzalini and Capitanio (1998)
# for notations and explanations

immutable MvSkewNormal <: Sampleable{Multivariate, Continuous}
    ξ::Vector{Float64}      # Location vector
    ω::Diagonal{Float64}      # Scale (diagonal matrix)
    Ωz::PDMat               # Correlation matrix
    α::Vector{Float64}      # Shape vector
    # Auxiliary parameters for sampling
    δ::Vector{Float64}
    Ωstar::PDMat  # For sampling...
    function MvSkewNormal(ξ::Vector{Float64}, Ω::Matrix{Float64}, α::Vector{Float64})
        ω, Ωz = scale_and_cor(PDMat(Ω))
        d = length(α)
        _δ = δ(Ωz, α)
        Ωstar = Array(Float64, d+1, d+1)
        Ωstar[:, 1] = [1; _δ]
        Ωstar[1,:] = [1 _δ']
        Ωstar[2:end, 2:end] = Ωz.mat
        Ωstar = PDMat(Ωstar)
        new(ξ, ω, Ωz, α, _δ, Ωstar)
    end
end

# Convenience constructor - no location parameters
MvSkewNormal(Ω::Matrix{Float64}, α::Vector{Float64}) = MvSkewNormal(zeros(length(α)), Ω, α)

function show(io::IO, dist::MvSkewNormal)
    println(io, "ξ = $(dist.ξ)")
    println(io, "Ω = ")
    println(io, dist.ω * dist.Ωz.mat * dist.ω)
    println(io, "α = $(dist.α)")
end

length(dist::MvSkewNormal) = length(dist.α)

function _rand!{T<:Real}(dist::MvSkewNormal, out::AbstractVector{T})
    xx = randn(length(out) + 1)
    unwhiten!(dist.Ωstar, xx)
    x0, x = xx[1], xx[2:end]
    copy!(out, dist.ω * (x0 > 0 ? x : -x) + dist.ξ)
end

mean(dist::MvSkewNormal) = sqrt(2/π)*dist.ω*dist.δ + dist.ξ
var(dist::MvSkewNormal) = (dist.ω^2)*(diag(dist.Ωz.mat) - (sqrt(2/π)*dist.δ).^2)

function cov(dist::MvSkewNormal)
    μz=sqrt(2/π)*dist.δ
    dist.ω * (dist.Ωz.mat - μz*μz') * dist.ω
end

# Special functions (related to cumulant function of half-normal distribution)
ζ₀(x::Float64) = log(2π * normcdf(x))
ζ₀(x::Vector{Float64}) = map(ζ₀, x)

ζ₁(x::Float64) = normpdf(x)/normcdf(x)
ζ₁(x::Vector{Float64}) = map(ζ₁, x)

# Note that these fit functions are not completely reliable.
# Sometimes a numerical errors occur, for instance, if the line
# search starts with too large a step size, then the function
# evaluations are undefined. On other occasions, the algorithm
# may converge only to a local optimum. The latter problem could
# perhaps be overcome through the use of the NLOpt library
# which allows the user to randomly generate multiple start points.
# The former problem could perhaps be resolved with error handling.

# Fit a MvSkewNormal regression model for model to data
# X = (n x p) model matrix (each row is corresponds to the predictors of one response)
# Y = (n x k) response matrix (each row correspond to one output observation)
#
# The following model is assumed:
#   yᵢ ∼ SNₖ(ξᵢ, Ω, α) where ξᵢ = xᵢ̱β
# for some (p x k) matrix β of parameters
function fit_MvSkewNormal(X::Matrix{Float64}, Y::Matrix{Float64}; kwargs...)
    size(X,1) == size(Y,1) || throw(ArgumentError("X and Y must have the same number of rows"))
    n,p = size(X)
    k = size(Y,2)
    # β is a (p x k) parameter matrix
    u(β::Matrix{Float64}) = Y - X*β
    V(β::Matrix{Float64}) = (u(β)'u(β))/n

    function ll(params::Vector{Float64})
        β = reshape(params[1:(p*k)], p, k)
        η = params[p*k+1:end]
        -(-0.5 * n * log(det(V(β))) - 0.5 * n * k + sum(ζ₀(u(β)*η)))
    end

    function dll!(params::Vector{Float64}, grad::Vector{Float64})
        β = reshape(params[1:(p*k)], p, k)
        η = params[p*k+1:end]

        uβ = u(β)
        Vβ = V(β)

        # ∂l∂β = X'u(β)*inv(Vβ) - X'ζ₁(u(β)*η)*η'

∂l∂β = X'*(Vβ\(uβ'))' - X'ζ₁(uβ*η)*η'

        ∂l∂η= uβ'ζ₁(u*η)
        grad[1:p*k] = -vec(∂l∂β)
        grad[p*k + 1:end] = -∂l∂η
    end

    function ll_and_dll!(params::Vector{Float64}, grad::Vector{Float64})
        β = reshape(params[1:p*k], p, k)
        η = params[p*k+1:end]

        println("β = $(β)")
        println("η = $(η)")

        uβ = u(β)
        Vβ = V(β)

        #
∂l∂β = X'u(β)*inv(V(β)) - X'ζ₁(u(β)*η)*η'

∂l∂β = X'*(Vβ\(uβ'))' - X'ζ₁(uβ*η)*η'
        ∂l∂η = uβ'ζ₁(uβ*η)
        grad[1:p*k] = -vec(∂l∂β)
        grad[p*k + 1:end] = -∂l∂η

        -(-0.5 * n * log(det(Vβ)) - 0.5 * n * k + sum(ζ₀(uβ*η)))
    end

    func = DifferentiableFunction(ll, dll!, ll_and_dll!)
    init = ones(p*k + k)

    results = optimize(func, init; kwargs...)
    β = reshape(results.minimum[1:p*k], p, k)
    η = results.minimum[p*k+1:end]
    Ω = V(β)
    ω = diagm(sqrt(diag(Ω)))
    α = ω * η

    return MvSkewNormal(Ω, α), β
end

# Fit a multivariate Skew Normal distribution
# directly to a set of observations
#
# Y = (n x k) observation matrix (each row correspond to one output observation)
function fit_MvSkewNormal(Y::Matrix{Float64}; kwargs...)
    dist, β = fit_skew(ones(size(Y,1),1), Y; kwargs...)
    return MvSkewNormal(vec(β), dist.ω*dist.Ωz.mat*dist.ω, dist.α)
end


text:
using Documenter
using TrajectoryOptimization
using RigidBodyDynamics

makedocs(
    sitename = "TrajectoryOptimization",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Introduction" => "index.md",
        "models.md",
        "costfunctions.md",
        "constraints.md",
        "problem.md",
        "solvers.md",
        "solving.md",
        "advanced.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/RoboticExplorationLab/TrajectoryOptimization.jl.git",
    devbranch = "master",
)


text:
module TestAmbiguities

include("preamble.jl")

using PkgWithAmbiguities

@testset begin
    @info "↓↓↓ Following failures are expected. ↓↓↓"
    results = @testtestset begin
        Aqua.test_ambiguities(PkgWithAmbiguities)
    end
    @info "↑↑↑ Above failures are expected. ↑↑↑"
    @test length(results) == 1
    @test results[1] isa Test.Fail

    # It works with other tests:
    Aqua.test_unbound_args(PkgWithAmbiguities)
    Aqua.test_undefined_exports(PkgWithAmbiguities)
end

end  # module


text:
module DO

using LinearAlgebra

export integrate, DOProblem

include("utils.jl")
include("low_rank_timestepping.jl")
mutable struct DOProblem
    flux # semidiscretization
    factorized # does the flux factorize? true/false
    U0 # initial set of modes n × s
    Z0 # initial set of stochastic coefficients l × s
end

function integrate(prob::DOProblem, ts, σ_min, integration_method = :direct_time_marching)
    U, Z = prob.U0, prob.Z0
    Us, Zs = [U], [Z]
    integrator = getfield(DO, Symbol(string(integration_method,:!)))
    for i in 2:length(ts)
        println(i)
        Δt = ts[i] - ts[i-1]
        if prob.factorized
            L_U, L_Z = prob.flux(U,Z,ts[i])
            N = Δt*normal_direction(U, Z, L_U, L_Z)
            if norm(N) > σ_min
                U, Z = augment_rank(N, U, Z, σ_min)
                println("   rank augmented")
            else
                integrator(Δt, U, Z, L_U, L_Z)
            end
        else
            L = prob.flux(U*Z', ts[i])
            N = Δt*normal_direction(U, Z, L)
            if norm(N) > σ_min
                U, Z = augment_rank(N, U, Z, σ_min)
                println("   rank augmented!")
            else
                integrator(Δt, U, Z, L)
            end
        end
        push!(Us, U)
        push!(Zs, Z)
    end
    return Us, Zs
end

end


text:
using DoubleFloats, GenericLinearAlgebra, LinearAlgebra
using Test

using Random
phi = Base.MathConstants.golden

include("matmul.jl")
include("specialvalues.jl")
include("arithmetic.jl")
include("functions.jl")

include("concrete_accuracy.jl")
include("function_accuracy.jl")

include("coverage.jl")
include("promotion.jl")

include("parse.jl")


text:
using Regression
using Base.Test

import Regression: loss, has_bias, pred_without_bias, pred_with_bias, initsol

d = 5
k = 3
n = 100
X = randn(d, n)
y = randn(n)
Y = randn(k, n)
c = rand(1:k, n)
b = 2.0

pb = linearreg(X, y; bias=b)
@test isa(pb, UnivariateRegression{SqrLoss})
@test has_bias(pb)
@test loss(pb) == SqrLoss()
@test pred_without_bias(pb) == LinearPred(d)
@test pred_with_bias(pb) == AffinePred(d, b)
@test initsol(pb) == zeros(d+1)

pb = linearreg(X, Y; bias=b)
@test isa(pb, MultivariateRegression{SumSqrLoss})
@test has_bias(pb)
@test loss(pb) == SumSqrLoss()
@test pred_without_bias(pb) == MvLinearPred(d, k)
@test pred_with_bias(pb) == MvAffinePred(d, k, b)
@test initsol(pb) == zeros(k, d+1)

pb = logisticreg(X, sign(y); bias=b)
@test isa(pb, UnivariateRegression{LogisticLoss})
@test has_bias(pb)
@test loss(pb) == LogisticLoss()
@test pred_without_bias(pb) == LinearPred(d)
@test pred_with_bias(pb) == AffinePred(d, b)
@test initsol(pb) == zeros(d+1)

pb = mlogisticreg(X, c, k; bias=b)
@test isa(pb, MultivariateRegression{MultiLogisticLoss})
@test has_bias(pb)
@test loss(pb) == MultiLogisticLoss()
@test pred_without_bias(pb) == MvLinearPred(d, k)
@test pred_with_bias(pb) == MvAffinePred(d, k, b)
@test initsol(pb) == zeros(k, d+1)


text:
## HierarchicalDomain

function Base.show(io::IO, H::HierarchicalDomain{S,T}) where {S,T}
    print(io,"$(nlevels(H))-level HierarchicalDomain{$S,$T}:\n")
    show(io,UnionDomain(collectdata(H)))
end

## HierarchicalSpace

function Base.show(io::IO, H::HierarchicalSpace{S,T}) where {S,T}
    print(io,"$(nlevels(H))-level HierarchicalSpace{$S,$T}:\n")
    show(io,PiecewiseSpace(collectdata(H)))
end

## HierarchicalMatrix{F<:GreensFun,G<:GreensFun}
# Base.writemime because HierarchicalMatrix <: AbstractArray

## HierarchicalOperator{U<:Operator,V<:AbstractLowRankOperator}

Base.alignment(io::IO, x::Infinity) = (1,0)
function Base.show(io::IO, ::MIME"text/plain", H::HierarchicalMatrix{F,GreensFun{L,T}}) where {F<:GreensFun,L<:LowRankFun,T}
    print(io,"$(nlevels(H))-level HierarchicalMatrix of GreensFun's with blockwise ranks:\n")
    Base.print_matrix(io,blockrank(H),"["," ","]")
end
function Base.show(io::IO, ::MIME"text/plain", H::HierarchicalOperator{U,V}) where {U<:Operator,V<:AbstractLowRankOperator}
    print(io,"$(nlevels(H))-level HierarchicalOperator with blockwise ranks:\n")
    Base.print_matrix(io,blockrank(H),"["," ","]")
end


text:
######################################################
# EXEMPLO SIMPLES DE USO DO PACOTE FLUX
# para treinamento de uma rede neural pequena para
# o problema de classificação de caracteres numéricos
# da MNIST.
#
# (adaptado do exemplo do "model-zoo" https://github.com/FluxML/model-zoo/blob/master/vision/mlp_mnist/mlp_mnist.jl)
#
# Autor: Leonardo D. Secchin
# Data : 31/08/2021
# leonardosecchin.github.io/topicospo/
#
# Execução:
#   treinar()
# ou
#   treinar(ninter = 30, eta = 1e-2, batchsize = 256, nepochs = 15)
######################################################

using Plots, Statistics

# carrega Flux
using Flux

# exporta estruturas internas do Flux para ambiente do usuário
# (assim não precisa escrever Flux.ESTRUTURA)
using Flux.Data: DataLoader
using Flux: onehotbatch, onecold, @epochs
using Flux.Losses: mse

# carrega pacote contendo base de dados MNIST
using MLDatasets

##########################
# CARREGA DATASET MNIST
##########################
function carrega_dados(batchsize)

    # carrega/baixa MNIST
    #   (xtrain,ytrain): pares entrada-saída para treino da rede
    #   (xtest,ytest)  : pares entrada-saída para teste da rede treinada
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest,  ytest  = MLDatasets.MNIST.testdata(Float32)

    # Reorganiza imagens para estrutura Flux
    xtrain = Flux.flatten(xtrain)
    xtest = Flux.flatten(xtest)

    # Codifica as saídas (labels '0' a '9') em vetores de 0's e 1's de tamanho 10 (one-hot-encode)
    #   Para visualizar o resultado, execute por exemplo  onehotbatch([1 2 3], 0:9)
    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    # Divide os dados de treinamento em mini-lotes de tamanho batchsize
    #   shuffle=true diz que os dados são embaralhados antes de cada época
    #   (veja pgs 8 e 13 de https://leonardosecchin.github.io/files/topicospo/ML_4-sg.pdf)
    #   DataLoader é uma estrutura própria do pacote Flux.jl. Execute ?DataLoader para detalhes.
    train_loader = DataLoader((xtrain, ytrain), batchsize=batchsize, shuffle=true)

    # Divide os dados de teste em mini-lotes
    test_loader = DataLoader((xtest, ytest), batchsize=batchsize)

    # retorna os dados de treino e teste
    return train_loader, test_loader
end

##########################
# MONTA A REDE NEURAL
##########################
function constroi_rede(ninter)
    # Chain e Dense são funções do pacote Flux para criação da rede neural
    #
    #   Dense(n[l-1],n[l], ativação) cria uma camada tradicional (como vista em aula)
    #   com n[l-1] neurônios entrando, n[l] saindo, e onde o fluxo é guiado por
    #   a[l] = a.(W[l] * a[l-1] .+ bias[l])
    #   Pode-se indicar diferentes funções de ativação: sigmoid, relu ...
    #   Veja ?Dense e a ajuda do Flux para mais detalhes
    #
    #   Chain agrega as camadas em um único objeto

    # Abaixo, uma rede com a seguinte estrutura é criada:
    #   camada de entrada com 784=28*28 neurônios
    #   ninter neurônios na segunda camada, todos com ativação sigmoid
    #   camada de saída com 10 neurônios, todos com ativação sigmoid
    #
    # Pesos e vieses são iniciados randomicamente
    return Chain(
 	        Dense(784, ninter, sigmoid),
            Dense(ninter, 10, sigmoid)
            )
end

##########################
# AVALIA RISCO EMPÍRICO (feedforward) (sobre todos os dados de treinamento)
#
# Algumas funções de perda implementadas no Flux:
#   mse (+ usada aula) : (ŷ .- y) .^ 2
#   msle               : log.(ŷ .+ ϵ) .- log.(y .+ ϵ)) .^ 2
#   crossentropy       : -sum(y .* log.(ŷ .+ ϵ))
#   binarycrossentropy : -y * log(ŷ + ϵ) - (1 - y) * log(1 - ŷ + ϵ)
#
# Legenda
#   y = a[L] é o rótulo do dado de entrada
#   ŷ é a resposta da rede
#
# Muitas outras funções estão disponíveis. Veja
# https://fluxml.ai/Flux.jl/stable/models/losses
#
# Obs: é possível implementar sua própria função de perda.
# Veja exemplos do model-zoo.
##########################
function loss(data_loader, rede)
    acc = 0
    ls = 0.0f0
    num = 0

    # percorre todos os dados em data_loader
    for (x, y) in data_loader
        # pares entrada-saída do treinamento
        x, y = cpu(x), cpu(y)

        # respostas da rede (feedforward)
        yy = rede(x)

        # soma perda dos dados calculados
        ls  += mse(yy, y, agg=sum)

        # calcula taxa de acerto (número de respostas corretas)
        # Aqui, "onecold" decodifica o dado do vetor de 0's e 1's para rótulo de '0' a '9'
        acc += sum(onecold(cpu(yy)) .== onecold(cpu(y)))

        # conta dados de treinamento calculados
        num += size(x, 2)
    end

    # retorna média das perdas (risco empírico) e percentual de acertos
    return ls/num, acc/num
end

##########################
# FUNÇÃO PRINCIPAL
##########################
function treinar(; ninter = 32, eta = 4e-1, batchsize = 256, nepochs = 20)

    # PARÂMETROS
    # ninter    : número de neurônios na camada intermediária
    # eta       : taxa de aprendizagem
    # batchsize : tamanho do mini-lote
    # nepochs   : número de épocas

    # carrega base de dados MNIST e cria mini-lotes
    train_loader, test_loader = carrega_dados(batchsize)

    # constrói rede
    # Obs: "cpu" refere-se à execução em CPU normal.
    # Flux suporta processamento por GPU, veja documentação/exemplos.
    rede = constroi_rede(ninter) |> cpu

    # variáveis da rede a serem otimizadas
    wb = Flux.params(rede)

    ##############################################
    ######### Algoritmo para treinamento #########
    ##############################################
    # Descent(η = 0.1)              : gradiente estocástico básico (SGD)
    # Momentum(η = 0.01, ρ = 0.9)   : SGD com momento
    # RMSProp(η = 0.001, ρ = 0.9)
    # ADAM(η = 0.001, β::Tuple = (0.9, 0.999))
    # ADAGrad(η = 0.1)
    # ADADelta(ρ = 0.9)
    #
    # *** Não citados em aula: ***
    # Nesterov(η = 0.001, ρ = 0.9)  : método acelerado de Nesterov
    # RADAM(η = 0.001, β::Tuple = (0.9, 0.999)): Rectified ADAM
    # AdaMax(η = 0.001, β::Tuple = (0.9, 0.999))
    # AMSGrad(η = 0.001, β::Tuple = (0.9, 0.999))
    # NADAM(η = 0.001, β::Tuple = (0.9, 0.999))
    # ADAMW(η = 0.001, β::Tuple = (0.9, 0.999), decay = 0)
    # OADAM(η = 0.0001, β::Tuple = (0.5, 0.9))
    # AdaBelief(η = 0.001, β::Tuple = (0.9, 0.999))
    #
    # Legenda
    #   η: taxa de aprendizagem fixa (SGD) ou "de referência" (veja Adagrad, pg 16 de https://leonardosecchin.github.io/files/topicospo/ML_4-sg.pdf)
    #   ρ: balanço do momento (β da pg 20 de https://leonardosecchin.github.io/files/topicospo/ML_4-sg.pdf)
    #   β: parâmetros de momento (não visto com detalhes em aula)
    # consulte https://fluxml.ai/Flux.jl/stable/training/optimisers/ para detalhes

    opt = Descent(eta)
#     opt = ADAM(eta)

    # medidas para montagem de gráficos
    train_loss_hist = []
    train_acc_hist  = []
    test_loss_hist  = []
    test_acc_hist   = []

    ###############################
    ######### Treinamento #########
    ###############################

    println("Treinando a rede neural...")

    # loop épocas
    for epoch in 1:nepochs

        # percorre os mini-lotes de treinamento (época)
        for (x, y) in train_loader

            # dados de treinamento do mini-lote
            x, y = cpu(x), cpu(y)

            # computa gradiente no ponto wb corrente
            # (backpropagation nos dados do mini-lote atual)
            g = gradient(() -> mse(rede(x), y), wb)

            # realiza uma iteração do método, atualizando pesos e vieses da rede
            Flux.Optimise.update!(opt, wb, g)

        end

        # calcula risco empírico sobre todos os dados de treinamento
        train_loss, train_acc = loss(train_loader, rede)

        # Calcula risco empírico nos dados de teste
        # (somente para impressão na tela).
        # Poderia ser usado para parar prematuramente (dados de validação).
        test_loss,  test_acc  = loss(test_loader, rede)

        # imprimi dados da época
        println()
        println("Época: $epoch")
        println("  Dados treinamento:  risco empir = $train_loss ,  % acertos = $train_acc")
        println("  Dados teste      :  risco empir = $test_loss ,  % acertos = $test_acc")

        # guarda medidas para montagem de gráficos
        append!(train_loss_hist, train_loss)
        append!(train_acc_hist, train_acc)
        append!(test_loss_hist, test_loss)
        append!(test_acc_hist, test_acc)
    end

    # Plota gráficos
    fig_loss = plot(title="Perda", xlabel="épocas", ylabel="", fontfamily="Computer Modern")
    fig_loss = plot!(1:length(train_loss_hist), train_loss_hist, lw=1, label="Dados treinamento (risco empírico)")
    fig_loss = plot!(1:length(test_loss_hist), test_loss_hist, lw=1, label="Dados teste")
    savefig(fig_loss, "fig_loss.pdf")

    fig_acc = plot(title="% acerto", xlabel="épocas", ylabel="% acerto", fontfamily="Computer Modern", legend=:bottomright)
    fig_acc = plot!(1:length(train_acc_hist), train_acc_hist, lw=1, label="Dados treinamento")
    fig_acc = plot!(1:length(test_acc_hist), test_acc_hist, lw=1, label="Dados teste")
    savefig(fig_acc, "fig_acc.pdf")
end

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
using Pixell
using Documenter

ENV["PLOTS_DEFAULT_BACKEND"] = "GR"
ENV["GKSwstype"] = "100"
using Plots
Plots.default(fontfamily="Computer Modern", fmt=:svg)

DocMeta.setdocmeta!(Pixell, :DocTestSetup,
    :(
      using Pixell, Plots);
    recursive=true)

makedocs(;
    modules=[Pixell],
    authors="Zack Li, Yilun Guan",
    repo="https://github.com/simonsobs/Pixell.jl/blob/{commit}{path}#{line}",
    sitename="Pixell.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://simonsobs.github.io/Pixell.jl",
        assets=["assets/custom.css"]
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Map Manipulation" => "map_manipulation.md",
        "Developer Notes" => "developer_notes.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/simonsobs/Pixell.jl",
    devbranch="main",
)


text:
# Unit tests for integrations/atomsbase_integration.jl

@testset "integrations/atomsbase_integration.jl" begin
    struct MockSystem1D <: AbstractSystem{1} end
    AtomsBase.bounding_box(s::MockSystem1D) = infinite_box(3)
    AtomsBase.boundary_conditions(s::MockSystem1D) = [Periodic(), Periodic(), Periodic()]
    Base.length(s::MockSystem1D) = 3
    AtomsBase.species_type(s::MockSystem1D) = Symbol
    Base.getindex(s::MockSystem1D, index::Int) = Symbol(index)
    AtomsBase.position(s::MockSystem1D) = [11, 12, 13]u"m"
    AtomsBase.atomic_symbol(s::MockSystem1D) = [Symbol(11), Symbol(12), Symbol(13)]
    AtomsBase.atomic_number(s::MockSystem1D) = [10, 10, 10]
    AtomsBase.atomic_mass(s::MockSystem1D) = [10, 10, 10]u"kg"
    AtomsBase.velocity(s::MockSystem1D) = missing
    AtomsBase.position(s::MockSystem1D, i) = i * u"m"
    AtomsBase.atomic_symbol(s::MockSystem1D, i) = Symbol(i)
    AtomsBase.atomic_number(s::MockSystem1D, i) = 0
    AtomsBase.atomic_mass(s::MockSystem1D, i) = 0u"kg"
    AtomsBase.velocity(s::MockSystem1D, i) = missing

    target = DynamicSystem(MockSystem1D(), 3.14u"s")

    @test bounding_box(target) == infinite_box(3)
    @test boundary_conditions(target) == [Periodic(), Periodic(), Periodic()]
    @test periodicity(target) == [true, true, true]
    @test n_dimensions(target) == 1
    @test length(target) == 3
    @test size(target) == (3,)
    @test species_type(target) == Symbol
    @test [x for x ∈ target] == [Symbol(1), Symbol(2), Symbol(3)]
    @test position(target) == [11, 12, 13]u"m"
    @test ismissing(velocity(target))
    @test AtomsBase.atomic_symbol(target) == [Symbol(11), Symbol(12), Symbol(13)]
    @test atomic_number(target) == [10, 10, 10]
    @test atomic_mass(target) == [10, 10, 10]u"kg"
    @test position(target, 1) == 1u"m"
    @test ismissing(velocity(target, 1))
    @test AtomsBase.atomic_symbol(target, 1) == Symbol(1)
    @test atomic_number(target, 1) == 0
    @test atomic_mass(target, 1) == 0u"kg"
end


text:
type SCM3GPPMacro
    # Angle spread
    mu_AS::Real
    eps_AS::Real
    r_AS::Real

    pathAS::Real

    # Delay spread
    mu_DS::Real
    eps_DS::Real
    r_DS::Real

    eps_PL::Real # 10*path loss exponent

    nPaths::Integer

    SCM3GPPMacro() = new()
end

function urbanMacro15Deg()
    scen = SCM3GPPMacro()

    # see Table 5.1 in 3GPP TR v12, p.17
    scen.mu_AS = 1.18
    scen.eps_AS = 0.21
    scen.r_AS = 1.3

    scen.pathAS = 2.0

    scen.mu_DS = -6.18
    scen.eps_DS = 0.18
    scen.r_DS = 1.7

    scen.eps_PL = 35
    scen.nPaths = 6

    return scen
end

function generate_channel(scen::SCM3GPPMacro, nAntennas; nCoherence=1,nBatches=1)
    h = zeros(Complex128,nAntennas, nCoherence, nBatches)
    t = zeros(Complex128,nAntennas, nBatches)
    for i in 1:nBatches

        # user angle
        theta = (rand()-0.5)*120

        # path delays
        DS = 10.^(scen.mu_DS + scen.eps_DS*randn());
        Tc = 1/3.84e6;
        tau = -scen.r_DS*DS*log.(rand(scen.nPaths));
        tau = sort(tau) - minimum(tau);
        tau_quant = Tc/16*floor.(tau/Tc*16 + 0.5);

        # path powers
        exponent = -1./DS*(scen.r_DS-1)/scen.r_DS;
        Z = randn(scen.nPaths)*3; #per path shadow fading in dB
        p = exp.(exponent*tau).*(10.^(0.1*Z));
        p = p./sum(p);

        # path AoDs
        AS = 10.^(scen.mu_AS + scen.eps_AS*randn());
        aodsm = randn(scen.nPaths)*scen.r_AS*AS;
        ixs = sortperm(abs.(aodsm));
        aodsm = aodsm[ixs];

        (h[:,:,i],t[:,i]) = scm_channel(theta .+ aodsm, p, nAntennas; nCoherence=nCoherence,AS=scen.pathAS)

        # pathloss
        minDist = 1000
        maxDist = 1500
        distance = rand() * (maxDist-minDist) + minDist
        PL = scen.eps_PL*log10(distance/maxDist)
        beta = 10^(-0.1*PL)

        h[:,:,i] .*= sqrt(beta)
        t[:,i] .*= beta
    end
    return (h,t)
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
# note: the functions (in cross_validation.jl have major redundancies with runmodel_multi.jl
# we may want to improve this

function dtm(dtmtable::DTMTable, settingsVector::Vector{ModelSettings};file::String=joinpath(mktempdir(), defaultModelNameWtihCSVext))
    # set rnd state to make this function reporducible (irrespective of trn/val idx)
    # Random.srand should not depend on sett.seed as we do not 'store' the original seed in the resulting Excel file.
    nSettings = length(settingsVector)

    intDatahash = floor(Int, .25 * hash(2231, hash(dtmtable.features, hash(dtmtable.numerator, hash(dtmtable.denominator, hash(dtmtable.weight))))))
    Random.seed!(intDatahash)

   # initialize variables
    path_and_fn_wo_extension = "some"
    path_and_fn_wo_extension, ext = splitext(file)
    header = Vector{String}(undef, 1)
    header_settings = Vector{String}(undef, 1)
    i = 1

        # run an initial dummy model to get the 'structure' (ie header and size) of the output statistics
    dummySett = deepcopy(settingsVector[1])
    dummySett.print_details = false;dummySett.writeTree = false;dummySett.writeSasCode = false;dummySett.writeIterationMatrix = false;dummySett.writeResult = false;dummySett.writeStatistics = false;dummySett.boolCreateZipFile = false;dummySett.writeCsharpCode = false;dummySett.writeVbaCode = false;dummySett.saveJLDFile = false;dummySett.saveResultAsJLDFile = false;dummySett.showProgressBar_time = false;dummySett.prroduceEstAndLeafMatrices = false;dummySett.write_dot_graph = false;
    path_and_fn_wo_extension_mod = string(path_and_fn_wo_extension, "_0")
    fnmod = string(path_and_fn_wo_extension_mod, ext)
    dummySett.iterations = min(10, dummySett.iterations) # is this a good idea? I think we should be fine.
    somestrings, model = run_model_actual(dtmtable, dummySett, fnmod)
    desc, numbrs, desc_settingsvec, settingsvec, selectedPerformanceMeasure = get_stats(model, perfMeasure=settingsVector[1].performanceMeasure)
            # add first column
    pushfirst!(numbrs, i)
    pushfirst!(settingsvec, string(i))
    pushfirst!(desc, "Number")
    pushfirst!(desc_settingsvec, "Number")
    header = deepcopy(desc)
    header_settings = deepcopy(desc_settingsvec)
    allstats = Array{Float64,2}(undef, length(numbrs), nSettings)
    allsettings = Array{Any,2}(undef, length(settingsvec), nSettings)
    allmodels = Vector{Any}(undef, nSettings)
        # define default stats if the model fails
    defaulted_stats = deepcopy(numbrs)
    defaulted_settings = deepcopy(settingsvec)
    emptyModel = deepcopy(EmtpyDTModel())
    fill!(defaulted_stats, 0.0)

        # create a dictionary which contains all data
    di = Dict("ext" => ext, "emptyModel" => emptyModel, "defaulted_settings" => defaulted_settings, "header_settings" => header_settings, "header" => header, "intDatahash" => intDatahash, "settingsVector" => deepcopy(settingsVector), "dtmtable" => dtmtable, "path_and_fn_wo_extension" => path_and_fn_wo_extension, "defaulted_stats" => defaulted_stats)
        # send data to all workers: without this command we will have a lot of overhead to send the data to the processes; but the data is constant and only needs to be sent ONCE!)
    sendto_module(DecisionTrees, Distributed.workers(), local_data_dict=deepcopy(di))
        # if Distributed.nprocs()==1
    local_data_dict = di # this is needed in case there are no workers
        # end
        # run all models in parallel
        # @warn("this is currently rather slow as the local_data_dict might be transferred to each worker for each iteration -> improve this!.... ? global const variable....?")
    pmapresult = Distributed.pmap(iLoop->singleRunDtm(iLoop, local_data_dict), 1:length(settingsVector))

        # 2. run models
    for pmap_entry in pmapresult# this_sample in cvsampler
        numbrs = pmap_entry[1]
        settingsvec = pmap_entry[2]
        i = pmap_entry[3]
        mdl = pmap_entry[4]
            # append stats
        allstats[:,i] .= deepcopy(numbrs)
        allsettings[:,i] .= deepcopy(settingsvec)
        allmodels[i] = mdl
    end

        # 3. aggregate some statistics
    statsdf = DataFrame(transpose(allstats), :auto)
    settsdf = DataFrame(permutedims(allsettings, [2,1]), :auto)
    DataFrames.rename!(settsdf, Symbol.(header_settings))

    fld, namestr = splitdir(path_and_fn_wo_extension)
    filen = string(path_and_fn_wo_extension, "_multistats.xlsx")
    if isfile(filen)
        @info "Deleting $(filen)."
        rm(filen)
    end

        # calc some 'stats of stats'
    col_rng_stats = 2:size(statsdf, 2)
    means = map(x->Statistics.mean(statsdf[!,x]), col_rng_stats)
    medians = map(x->Statistics.median(statsdf[!,x]), col_rng_stats)
    mins = map(x->minimum(statsdf[!,x]), col_rng_stats)
    maxs = map(x->maximum(statsdf[!,x]), col_rng_stats)
    stds = map(x->StatsBase.std(statsdf[!,x]), col_rng_stats)
    sums = map(x->sum(statsdf[!,x]), col_rng_stats)

    stats_of_stats = hcat(means, medians, mins, maxs, stds, sums)
    titles = ["Mean","Median","Min","Max","Std Dev","Sum"]
    stats_of_stats = hcat(titles, transpose(stats_of_stats))
    stats_of_stats = vcat(repeat([""], 1, size(stats_of_stats, 2)), stats_of_stats)

    allstats_with_stats = vcat(transpose(allstats), stats_of_stats)
        # NOTE! statsdf is RE defined here!
    statsdf = DataFrame(allstats_with_stats, :auto)
    DataFrames.rename!(statsdf, Symbol.(header))

        # define Exceldata
    sh1 = ExcelSheet("settings", settsdf)
    sh2 = ExcelSheet("stats", statsdf)
    xld = ExcelData(ExcelSheet[sh1,sh2], Array{Chart,1}(undef, 0))
        # write data
    try
        println("Writing summary statistics to file:\r\n$(filen)")
        @time writeStatistics(xld, filen, true, false)
    catch e
        @show e
        @warn("DTM: Failed to create Excel Statistics file. You may want to check the PyCall installation and whether the required Python packages are installed. \r\n $(filen)")
    end

    return statsdf, settsdf, allmodels
end

function singleRunDtm(i::Int, local_data_dict::Dict)
    di = local_data_dict
    defaulted_settings = di["defaulted_settings"]
    defaulted_stats = di["defaulted_stats"]
    emptyModel = di["emptyModel"]

    try
        ext = di["ext"]
        header_settings = di["header_settings"]
        header = di["header"]
        intDatahash = di["intDatahash"]
        settingsVector = di["settingsVector"]
        dtmtable = di["dtmtable"]
        path_and_fn_wo_extension = di["path_and_fn_wo_extension"]

        sett = settingsVector[i]

        # set seed
        sett.seed = floor(Int, 0.31 * hash(intDatahash, hash(99812, hash(i))))
        # define path
        path_and_fn_wo_extension_mod = string(path_and_fn_wo_extension, "_", i)
        fnmod = string(path_and_fn_wo_extension_mod, ext)
        # run model
        somestrings, model = run_model_actual(dtmtable, sett, fnmod)
        desc, numbrs, desc_settingsvec, settingsvec, selectedPerformanceMeasure = get_stats(model, perfMeasure=sett.performanceMeasure)
        # add index as first column
        pushfirst!(numbrs, i)
        pushfirst!(settingsvec, string(i))
        pushfirst!(desc, "Number")
        pushfirst!(desc_settingsvec, "Number")
        # check if header and settings header are as expected
        if !all(header[2:end] .== desc[2:end])
            @warn("Model run $(i) returned an unexpected results vector:")
            @show desc
            @show header
        end
        if !all(header_settings[2:end] .== desc_settingsvec[2:end])
            @warn("Model run $(i) returned an unexpected results vector:")
            @show desc_settingsvec
            @show header_settings
        end
        return numbrs, settingsvec, i, model
    catch eri
        @warn("DTM: Model $(i) failed.")
        println(eri)
        return defaulted_stats, defaulted_settings, i, emptyModel
    end

    return nothing
end	 # singleRunDtm

function get_model_stats(resulting_model::DTModel, defaulted_modelstats_df::DataFrame, i::Int)
	thesemodelstats = deepcopy(resulting_model.modelstats)
	nm = DataFrame(ones(Int, size(thesemodelstats, 1), 1) .* i, :auto)
	DataFrames.names!(nm, Symbol[symbol("ModelNumber")])
	return  thesemodelstats = hcat(nm, thesemodelstats)
end

function createGridSearchSettings(sett::ModelSettings;args...)
    settingsVector = [sett]
	seen = []
	for (smember, valuelist) in args
		if in(smember, seen)
			error("DTM: (createGridSearchSettings) You provided multiple values for the field $(smember)")
		end
        push!(seen, smember)
        try
            sz = size(valuelist)
            @assert length(sz) == 1 "DTM: expected a vector for $(smember). You provided $(typeof(valuelist)). \r\n$(valuelist)"
            n = sz[1]
            @assert n > 0

            settingsVectorTemplate = deepcopy(settingsVector)

            for i = 1:n
                # set value for this 'batch'
                v = valuelist[i]
                thisSettingsBatch = deepcopy(settingsVectorTemplate)
                for kk = 1:length(thisSettingsBatch)
                    s = thisSettingsBatch[kk]
                    # updateSettingsMod!(s,smember=)
                    oldValue = getfield(s, smember)
                    if typeof(v) != typeof(oldValue)
                        valConverted = convertFromString(oldValue, string(v)) # tbd maybe string(v) could/should be replaced with v here
                        setfield!(s, smember, valConverted)
                    else
                        setfield!(s, smember, v)
                    end
                    res = checkIfSettingsAreValid(s)
                    @assert res == nothing "DTM: Some settings are invalid! Abort."
                end
                if i == 1
                    settingsVector = deepcopy(thisSettingsBatch)
                else
                    append!(settingsVector, thisSettingsBatch)
                end
            end
        catch e
            @show smember, size(valuelist)
            @show valuelist
            @show e
            throw("DTM: unable to create settings vector.")
        end
	end
	@assert length(unique(map(x->x.performanceMeasure, settingsVector))) == 1 "DTM: Perf measure must be the same for the entire settings Vector"
	println("DTM: Size of settings vector is $(length(settingsVector))")
    return settingsVector
end


text:
# -------------------------------------------------------------------------------
# Projection PatMat hinge
# -------------------------------------------------------------------------------
# general solver
function projection_exact(α0::AbstractVector,
                          β0::AbstractVector,
                          δ0::Real,
                          C1::Real,
                          C2::Real;
                          solver::Any = Optimizer(verbose = false))

    α = Convex.Variable(length(α0))
    β = Convex.Variable(length(β0))
    δ = Convex.Variable()

    objective   = Convex.sumsquares(α - α0)/2 + Convex.sumsquares(β - β0)/2 +
                  Convex.sumsquares(δ - δ0)/2
    constraints = [Convex.sum(α) == Convex.sum(β),
                   α <= C1,
                   α >= 0,
                   β <= C2*δ,
                   β >= 0]

    problem = Convex.minimize(objective, constraints)
    Convex.solve!(problem, solver)

    return vec(α.value), vec(β.value), δ.value
end

# our solution
function f(μ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real)
    return C2*δ0 + C2^2*sum(max.(β0 .- μ, 0)) - μ
end

function h(μ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real)
    λ = f(μ, α0, β0, δ0, C1, C2)
    return sum(min.(max.(α0 .- λ, 0), C1)) - sum(min.(max.(β0 .+ λ, 0), λ + μ))
end

function projection(α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real)

    if δ0 <= - C2*sum(max.(β0 .+ maximum(α0), 0))
        return zero(α0), zero(β0), zero(δ0)
    end

    μ_ub1 = minimum(β0)
    μ_ub2 = minimum(C2*δ0 + C2^2*sum(β0) .- α0)/(C2^2*length(β0) + 1)
    μ_ub  = min(μ_ub1, μ_ub2) - 1e8
    μ_lb  = maximum(β0) + C2*max(δ0, 0) + 1e8

    μ = Roots.find_zero(μ -> h(μ, α0, β0, δ0, C1, C2), (μ_lb, μ_ub))
    λ = f(μ, α0, β0, δ0, C1, C2)

    α = @. min(max(α0 - λ, 0), C1)
    β = @. min(max(β0 + λ, 0), λ + μ)
    δ = (λ + μ)/C2

    return α, β, δ
end

# -------------------------------------------------------------------------------
# Projection PatMat quadratic
# -------------------------------------------------------------------------------
# general solver
function projection_exact(α0::AbstractVector,
                          β0::AbstractVector,
                          δ0::Real;
                          solver::Any = Optimizer(verbose = false))

    α = Convex.Variable(length(α0))
    β = Convex.Variable(length(β0))
    δ = Convex.Variable()

    objective   = Convex.sumsquares(α - α0)/2 + Convex.sumsquares(β - β0)/2 +
                  Convex.sumsquares(δ - δ0)/2
    constraints = [Convex.sum(α) == Convex.sum(β),
                   α >= 0,
                   β >= 0,
                   δ >= 0]

    problem = Convex.minimize(objective, constraints)
    Convex.solve!(problem, solver)

    return vec(α.value), vec(β.value), δ.value
end

# our solution
function f(λ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real)
    return sum(max.(α0 .- λ, 0)) - sum(max.(β0 .+ λ, 0))
end

function projection(α0::AbstractVector, β0::AbstractVector, δ0::Real)

    if - maximum(β0) > maximum(α0)
        return zero(α0), zero(β0), zero(δ0)
    end

    λ = Roots.find_zero(λ -> f(λ, α0, β0, δ0), (- maximum(β0), maximum(α0)))

    α = @. max(α0 - λ, 0)
    β = @. max(β0 + λ, 0)
    δ = max(δ0, 1e-4)

    return α, β, δ
end

# -------------------------------------------------------------------------------
# Projection TopPushK hinge
# -------------------------------------------------------------------------------
# general solver
function projection_exact(α0::AbstractVector,
                          β0::AbstractVector,
                          C::Real,
                          K::Integer;
                          solver::Any = Optimizer(verbose = false))

    α = Convex.Variable(length(α0))
    β = Convex.Variable(length(β0))

    objective   = Convex.sumsquares(α - α0)/2 + Convex.sumsquares(β - β0)/2
    constraints = [Convex.sum(α) == Convex.sum(β),
                   α <= C,
                   α >= 0,
                   β <= Convex.sum(α)/K,
                   β >= 0]

    problem = Convex.minimize(objective, constraints)
    Convex.solve!(problem, solver)

    return vec(α.value), vec(β.value)
end

# our solution
function f(μ::Real, s::AbstractVector, K::Integer)
    i, j, d  = 2, 1, 1
    λ, λ_old = s[1], 0
    g, g_old = -K*μ, -K*μ

    while g < 0
        g_old = g
        λ_old = λ

        if s[i] <= s[j] + μ
            g += d*(s[i] - λ)
            λ  = s[i]
            d += 1
            i += 1
        else
            g += d*(s[j] + μ - λ)
            λ  = s[j] + μ
            d -= 1
            j += 1
        end
    end
    return -(λ - λ_old)/(g - g_old)*g_old + λ_old
end

function h(μ::Real, s::AbstractVector, α0::AbstractVector, β0::AbstractVector, C::Real, K::Integer)
    λ = f(μ, s, K)
    return sum(min.(max.(α0 .- λ .+ sum(max.(β0 .+ λ .- μ, 0))/K, 0), C)) - K*μ
end

function projection(α0::AbstractVector, β0::AbstractVector, C::Real, K::Real)

    if K >= length(β0)
        @error "No feasible solution: K < length(β0) needed."
        return α0, β0
    end

    if mean(partialsort(β0, 1:K; rev = true)) + maximum(α0) <= 0
        return zero(α0), zero(β0)
    end

    s    = vcat(.- sort(β0; rev = true), Inf)
    μ_lb = 1e-10
    μ_ub = length(α0)*C/K + 1e-6

    μ = Roots.find_zero(μ -> h(μ, s, α0, β0, C, K), (μ_lb, μ_ub))
    λ = f(μ, s, K)
    δ = sum(max.(β0 .+ λ .- μ, 0))/K

    α = @. max(min(α0 - λ + δ, C), 0)
    β = @. max(min(β0 + λ, μ), 0)

    return α, β
end

# ---------------------------------------------------------------------------------------
# Projection TopPushK quadratic
# ---------------------------------------------------------------------------------------
# general solver
function projection_exact(α0::AbstractVector,
                          β0::AbstractVector,
                          K::Integer;
                          solver::Any = Optimizer(verbose = false))

    α = Convex.Variable(length(α0))
    β = Convex.Variable(length(β0))

    objective   = Convex.sumsquares(α - α0)/2 + Convex.sumsquares(β - β0)/2
    constraints = [Convex.sum(α) == Convex.sum(β),
                   α >= 0,
                   β <= Convex.sum(α)/K,
                   β >= 0]

    problem = Convex.minimize(objective, constraints)
    Convex.solve!(problem, solver)

    return vec(α.value), vec(β.value)
end

# our solution
function h(μ::Real, s::AbstractVector, α0::AbstractVector, β0::AbstractVector, K::Integer)
    λ = f(μ, s, K)
    return sum(max.(α0 .- λ .+ sum(max.(β0 .+ λ .- μ, 0))/K, 0)) - K*μ
end

function projection(α0::AbstractVector, β0::AbstractVector, K::Integer)

    if K >= length(β0)
        @error "No feasible solution: K < length(β0) needed."
        return α0, β0
    end

    if mean(partialsort(β0, 1:K; rev = true)) + maximum(α0) <= 0
        return zero(α0), zero(β0)
    end

    s    = vcat(.- sort(β0; rev = true), Inf)
    μ_lb = 1e-10
    μ_ub = length(α0)*(maximum(α0) + maximum(β0))/K

    μ = Roots.find_zero(μ -> h(μ, s, α0, β0, K), (μ_lb, μ_ub))
    λ = f(μ, s, K)
    δ = sum(max.(β0 .+ λ .- μ, 0))/K

    α = @. max(α0 - λ + δ, 0)
    β = @. max(min(β0 + λ, μ), 0)

    return α, β
end


text:
@testset "test_matrix_comps" begin
A = [-0.21 0.2; 0.2 -0.21]
B = 0.01*[1 0; 0 1]
C = [1 0; 0 1]
D = 0
sys = ss(A,B,C,D)
sysr, G = balreal(sys)

@test gram(sysr, :c) ≈ G
@test gram(sysr, :o) ≈ G
@test sort(pole(sysr)) ≈ sort(pole(sys))

sysb,T = ControlSystems.balance_statespace(sys)
Ab,Bb,Cb,T = ControlSystems.balance_statespace(A,B,C)

@test Ab*T ≈ T*A
@test Bb ≈ T*B
@test Cb*T ≈ C

@test sysb.A ≈ Ab

@test ControlSystems.balance_transform(A,B,C) ≈ ControlSystems.balance_transform(sys)

W = [1 0; 0 1]
@test covar(sys, W) ≈ [0.002560975609756 0.002439024390244; 0.002439024390244 0.002560975609756]
D2 = [1 0; 0 1]
@test covar(ss(A,B,C,D2, 1), W) ≈ [1.000110108378310 -0.000010098377310; -0.000010098377310 1.000110108378310]
# Direct term means infinite covariance
@test covar(ss(A,B,C,D2), W) ≈ [Inf Inf; Inf Inf]

# No noise on second output should give finite variance
@test covar(ss(A,B,C,[1 0; 0 0]), W) ≈ [Inf Inf; Inf 0.002560975609756]

# Unstable system has inf covar
@test covar(ss([1 0; 0 1],B,C,0), W) == [Inf Inf; Inf Inf]

# Discrete system can have direct term
@test covar(ss(A,B,C,D2,0.1),W) ≈ [1.00011010837831 -1.0098377309782909e-5; -1.0098377309782909e-5 1.00011010837831]

# TODO test in Julia 0.7 to see if supported
# # Test special matrices
As = sparse(A)
Bs = sparse(B)
Cs = sparse(C)
Asb,Bsb,Csb,Ts = ControlSystems.balance_statespace(As,Bs,Cs) #Error no LAPACK function
#
# @test Abs*Ts ≈ Ts*As
# @test Bbs ≈ Ts*Bs
# @test Cbs*Ts ≈ Cs

# Test special values
Ar = rationalize.(A)
Br = rationalize.(B)
Cr = rationalize.(Float64.(C))    # When did ever rationalize work on Int?
Arb,Brb,Crb,Tr = ControlSystems.balance_statespace(Ar,Br,Cr)

@test Arb*Tr ≈ Tr*Ar
@test Brb ≈ Tr*Br
@test Crb*Tr ≈ Cr

end


text:
x = 2_000_000
ba = BitArray(repeat([1], x))
ba[1] = 0

for i = 4:2:x
  ba[i] = 0
end

for i = 3:2:x
  for j = 2:Int(floor(x/i))
    ba[i*j] = 0
  end
end

t = 0
for i = 2:x
  global t
  if ba[i]
    t += i
  end
end

println(t)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
abstract type StepController end
struct ConstantSteps <: StepController end
function propose_step!(::ConstantSteps, integ)
    return integ.dt
end

struct StandardSteps <: StepController end
function propose_step!(::StandardSteps, integ)

    @unpack dt, EEst = integ
    @unpack q = integ.constants
    @unpack gamma, qmin, qmax = integ.opts

    if iszero(EEst)
        q = inv(qmax)
    else
        localconvrate = q+1
        qtmp = DiffEqBase.fastpow(EEst, 1/localconvrate) / gamma
        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
    end
    dt /= q
    return dt
end

struct PISteps <: StepController end
function propose_step!(::PISteps, integ)
    # PI-controller
    @unpack dt, EEst, qold = integ
    @unpack gamma, beta1, beta2, qmin, qmax = integ.opts

    if iszero(EEst)
        q = inv(qmax)
    elseif EEst < 1  # Accepted
        q11 = DiffEqBase.value(DiffEqBase.fastpow(EEst, beta1))
        q = q11 / DiffEqBase.fastpow(qold, beta2) / gamma

        @fastmath q = DiffEqBase.value(max(inv(qmax),min(inv(qmin), q)))
        integ.qold = max(integ.EEst, integ.opts.qoldinit)
    else # Rejected step
        q11 = DiffEqBase.value(DiffEqBase.fastpow(EEst, beta1))
        q = q11 / gamma
        @fastmath q = DiffEqBase.value(max(inv(qmax),min(inv(qmin), q)))
    end
    dt /= q
    return dt
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
using Weave

docs_path = splitpath(@__DIR__)
while docs_path[end] ≠ "SpinModels"
    pop!(docs_path)
end
iszero(length(docs_path)) && error("Cannot locate `SpinModels` dir!")
docs_path = joinpath(push!(docs_path, "docs"))

jmddir(x...) = joinpath(@__DIR__, x...)

weave(jmddir("energy.jmd"), out_path = docs_path, doctype = "github");
weave(jmddir("ratiomatch.jmd"), out_path = docs_path, doctype = "github");
weave(jmddir("samplers.jmd"), out_path = docs_path, doctype = "github");


text:
macro define(name, definition)
    quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

_big(x) = x
_big(x::Number) = big(x)
_big(x::String) = parse(BigFloat, x)

function _big(x::Expr)
    y = x
    y.args .= _big.(y.args)
    return  y
end

macro big(x)
    return esc(_big(x))
end


text:
using Documenter, AIBECS
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# Generate examples
include("generate.jl")

EXAMPLES_jl = [f for f in readdir(EXAMPLEDIR) if endswith(f, ".jl")]
EXAMPLES_md = [replace(f, ".jl" => ".md") for f in EXAMPLES_jl]
GENERATEDEXAMPLES = [joinpath("examples", "generated", f) for f in EXAMPLES_md]

makedocs(
    sitename="AIBECS.jl",
    doctest = false, # TODO guessing I should remove that when actually deploying?
    # options
    modules = [AIBECS],
    # organisation
    pages = Any[
        "Home" => "index.md",
        "Prerequisites" => "prerequisites.md",
        "Examples" => GENERATEDEXAMPLES,
        "Function index" => "functions.md"
        ]
)

# Deploy
deploydocs(
    repo = "github.com/briochemc/AIBECS.jl.git",
)