text:
############################# Functions used in the computation of the Transition (in addition to the ones in functions.jl) [measure_trans, dynamicproblem_trans_vec_t, static problem period 2 (alt timing), simulate_trans, transition_vec2 (NLSolve), transition_vec2 (standard)] ########################

########################## Measure Trans #############################
function KLS3_measure_trans(p,rt,N,M0)
@unpack z_grid,a_grid,z_P=p #allows to use the parameters in p without having to addend "p.". Furthermore, it makes clear which parameters are used in each function.

    numZ = length(z_grid)
    numA = length(a_grid)

    M = zeros(numA,numZ,N)
    M[:,:,1] = M0
    M[:,:,2] = M0 # in the second period news arrives after states are determined
    Mnew = copy(M0)

    P=z_P
    for n=2:N-1
        apt_ind = rt[n].ap_ind
        Mold = copy(Mnew)
        Mnew = zeros(numA,numZ)
        for j=1:numZ #Old productivity state
            PP=@view P[j:j,:] #@view (and @views) lets P[j:j,:] reference the variable P instead of allocating a new array, P[j:j,:]. This speeds up the code and is thus used several times.
            for i=1:numA #Old asset state
                @views Mnew[apt_ind[i,j]:apt_ind[i,j],:] += Mold[i,j]*PP
            end
        end

        M[:,:,n+1] = copy(Mnew)

    end
return M

end

####################### Static Problem Period 2 #####################

function KLS3_staticproblem_period2_altTiming(p_o_t,p)

@unpack a_grid,z_grid_size,z,α,α_m,τ,Yf,τ_x,σ,δ,F_base,a_grid_size,r = p #See comment at the beginning of KL3_measure_trans
    ### Useful objects

    #Assets
    p_o_t=merge(p_o_t,(a_grid_mat = a_grid'*ones(1,z_grid_size),

    #Productivity
    z_grid_mat = z))

    # capital gain
    cap_gain = (1-p_o_t.Pk/p_o_t.Pk_lag)

    # Capital
    p_o_t=merge(p_o_t,(k_x = rt[1].k,
    k_nx = rt[1].k))

    # Useful constants
    MC = (p_o_t.Pk/α_m )^α_m * (p_o_t.w /((1-α)*(1-α_m)))^((1-α)*(1-α_m))

### Exporters

    # useful constant
    const_μ_x = (1 ./p_o_t.z_grid_mat).^(1/(α.*(1-α_m))) .* 1 ./p_o_t.k_x .* ((σ-1)./(σ)).^σ .*(p_o_t.ϕ_h + (p_o_t.ξ/τ).^σ.*Yf.*τ.*(1+τ_x).^(-σ))

    # marginal cost
    p_o_t=merge(p_o_t,(μ_x = const_μ_x.^( (α.*(1-α_m)) ./ (1+(σ-1).*α.*(1-α_m)) ) .*MC.^(1 ./(1+(σ-1).*α.*(1-α_m))),))

    # variable inputs
    p_o_t=merge(p_o_t,(m_x = (α_m./p_o_t.Pk) .* (p_o_t.z_grid_mat.*p_o_t.μ_x).^(1 ./(α.*(1-α_m))) .* p_o_t.k_x .*MC.^(-1 ./(α.*(1-α_m))),
    n_x = ( ((1-α).*(1-α_m)) ./ p_o_t.w ) .* ( p_o_t.z_grid_mat.*p_o_t.μ_x ).^(1 ./(α.*(1-α_m))) .* p_o_t.k_x .*MC.^(-1 ./(α.*(1-α_m))),

    # output
    yd_x = ((σ-1)/σ).^σ .* p_o_t.ϕ_h .* p_o_t.μ_x.^(-σ),
    yf_x = ((σ-1)/σ).^σ .* p_o_t.ξ.^σ .* Yf .* (τ .*(1+τ_x) .*p_o_t.μ_x).^(-σ)))
    p_o_t=merge(p_o_t,(y_x = p_o_t.yd_x+τ*p_o_t.yf_x,

#Prices adjust (change in demand)
    pd_x = (p_o_t.yd_x./p_o_t.ϕ_h).^(-1/σ)))
    p_o_t.pd_x[p_o_t.yd_x.==0] .= 0
    p_o_t=merge(p_o_t,(pf_x = ((p_o_t.yf_x/Yf).^(-1/σ))/(1+τ_x),))
    p_o_t.pf_x[p_o_t.yf_x.==0] .= 0

#Profits (in units of the final good)
    p_o_t=merge(p_o_t,(π_x = p_o_t.pd_x.*p_o_t.yd_x .+ p_o_t.ξ*p_o_t.yf_x.*p_o_t.pf_x .- ((r+δ).+(1-δ)*cap_gain)*p_o_t.Pk_lag*p_o_t.k_x .- p_o_t.w*p_o_t.n_x .- p_o_t.Pk*p_o_t.m_x .- p_o_t.w*F_base,))

# Non-Exporters

    # useful constant
    const_μ_nx = (1 ./p_o_t.z_grid_mat).^(1/(α.*(1-α_m))) .* 1 ./p_o_t.k_nx .* ((σ-1)./(σ)).^σ .* p_o_t.ϕ_h

    # marginal cost
    p_o_t=merge(p_o_t,(μ_nx = const_μ_nx.^((α.*(1-α_m))./(1+(σ-1).*α.*(1-α_m)))*MC.^(1 ./(1 + (σ-1).*α.*(1-α_m))),))

    # variable inputs
    p_o_t=merge(p_o_t,(m_nx = (α_m/p_o_t.Pk) .* (p_o_t.z_grid_mat.*p_o_t.μ_nx).^(1 ./(α*(1-α_m))) .* p_o_t.k_nx .*MC.^(-1 ./(α.*(1-α_m))),
    n_nx = ( ((1-α).*(1-α_m)) ./ p_o_t.w ) .* ( p_o_t.z_grid_mat.*p_o_t.μ_nx ).^(1 ./(α.*(1-α_m))) .* p_o_t.k_nx .*MC.^(-1 ./(α.*(1-α_m))),

    # output
    yd_nx = ((σ-1)/σ).^σ .* p_o_t.ϕ_h .* p_o_t.μ_nx.^(-σ)))
    p_o_t=merge(p_o_t,(yf_nx = zeros(size(p_o_t.yd_nx)),))
    p_o_t=merge(p_o_t,(y_nx = (p_o_t.yd_nx+τ*p_o_t.yf_nx),

#Prices
    pd_nx = ((p_o_t.yd_nx/p_o_t.ϕ_h).^(-1/σ))))
	p_o_t.pd_nx[p_o_t.yd_nx.==0] .= 0
    p_o_t=merge(p_o_t,(pf_nx = zeros(a_grid_size,z_grid_size),

#Profits
    π_nx =  p_o_t.pd_nx.*p_o_t.yd_nx  .-  ((r+δ).+(1-δ)*cap_gain)*p_o_t.k_nx*p_o_t.Pk_lag .- p_o_t.w*p_o_t.n_nx .- p_o_t.Pk*p_o_t.m_nx,

# Export decision
    e=BitArray(undef,a_grid_size,z_grid_size))) #BitArray: Array of only 1's and 0's, uses less memory than standard Float64 arrays (hence, makes the code faster)
    p_o_t.e .= p_o_t.π_x .>= p_o_t.π_nx

    return p_o_t
end

############## Dynamic Problem Transition Vectorized t ###############

function KLS3_dynamicproblem_trans_vec_t(p_o_t,p,rt)

 @unpack N,a_grid_size,z_grid_size,a_grid,r,γ,z_P,a_grid_vec,β = p #See comment at the beginning of KL3_measure_trans

### Dynamic problem: value functions and policy functions

### Shocks

#Generate solution objects
    v_p=copy(rt[N].v)
    if s.transition_AD == 1 # See comment in line 202 of main.jl.
    v_new = zeros(eltype(p_o_t.w),size(v_p))
    else
    v_new = zeros(size(v_p))
    end

    ap_ind = zeros(Int64,size(v_new))
    ap=zeros(a_grid_size,z_grid_size)
    asset_income = a_grid'.*(1+r)
    exponentg=1-γ
### Dynamic problem: value functions and policy functions
    P=z_P
    v_pt=v_p'
    ones_vec=β*ones(a_grid_size,1)
    neg_c_indexes=BitArray(undef,a_grid_size,z_grid_size) #BitArray: Array of only 1's and 0's, uses less memory than standard Float64 arrays (hence, makes the code faster)

    #Value function iteration algorithm

    for t=N-1:-1:2
        profits = rt[t].profits
        for j = 1:z_grid_size

                c =  profits[:,j] + asset_income .- a_grid_vec
                neg_c_indexes = real(c).<=0 #indices of a' for which consumption<0
                u = (c.^exponentg)/exponentg + neg_c_indexes*-1e50
                MM=ones_vec*P[j,:]'
                v,index_ap = findmax(real(u + MM*v_pt),dims=2)
                v_new[:,j] = v
                for (i,index) in enumerate(index_ap)
                ap_ind[i,j] = index[2]
                end
                ap[:,j]=a_grid[ap_ind[:,j]]
        end

        v_pt=copy(v_new') #Store output from value function iteration
        rt[t]=merge(rt[t],(v = copy(v_new),
        ap = copy(ap),
        ap_ind= copy(ap_ind)))
        rt[t]=merge(rt[t],(c = profits + rt[t].a_grid_mat.*(1+r) .-rt[t].ap,

    ## Store output to be used in simulation

        pd = (1 .-rt[t].e).*rt[t].pd_nx .+ rt[t].e.*rt[t].pd_x,
        yd = (1 .-rt[t].e).*rt[t].yd_nx .+ rt[t].e.*rt[t].yd_x,
        pf = rt[t].e.*rt[t].pf_x,
        yf = rt[t].e.*rt[t].yf_x,
        k = (1 .-rt[t].e).*rt[t].k_nx .+ rt[t].e.*rt[t].k_x,
        n_x = rt[t].e.*(rt[t].n_x),
        n_nx = (1 .-rt[t].e).*rt[t].n_nx))
        rt[t]=merge(rt[t],(n = (1 .-rt[t].e).*rt[t].n_nx .+ rt[t].e.*rt[t].n_x,
        m = (1 .-rt[t].e).*rt[t].m_nx .+ rt[t].e.*rt[t].m_x))

        rt[t]=merge(rt[t],(pd_nx = nothing,
        pd_x = nothing,
        yd_nx=nothing,
        yd_x=nothing,
        pf_x=nothing,
        yf_x=nothing,
        k_nx =nothing,
        k_x = nothing,
        n_nx=nothing,
        n_x=nothing,
        m_nx=nothing,
        m_x=nothing))

        rt[t]=merge(rt[t],(π_1 = (1 .-rt[t].e).*rt[t].π_nx .+ rt[t].e.*rt[t].π_x,
        a = copy(rt[t].a_grid_mat)))

    end

return rt
end
############## Simulate Transition ###############

function KLS3_simulate_trans(p,p_o_t,rt,Guess)
@unpack N,ω_h_c,σ,ω_h_k,rv,Pfv,Yfv,θ_v,β_v,δ_v,τ_v,τ_x_v,τ_m_c_v,τ_m_k_v,δ,Pm_k,τ_m_k,ω_m_k,Pm_c,τ_m_c,ω_m_c,r = p #See comment at the beginning of KL3_measure_trans

    wguess = [p_o_t.wt[1] exp.(Guess[N-1:2*N-4])' p_o_t.wt[N]]
    ξguess = [p_o_t.ξt[1] exp.(Guess[2*N-3:3*N-6])' p_o_t.ξt[N]]
    Pkguess = [p_o_t.Pkt[1] exp.(Guess[3*N-5:4*N-8])' p_o_t.Pkt[N]]

    if s.tariffsincome == 1
        Ycguess = [p_o_t.Yct[1] exp.(Guess[1:N-2])' p_o_t.Yct[N]]
        Ykguess = [p_o_t.Ykt[1] exp.(Guess[4*N-7:5*N-10])' p_o_t.Ykt[N]]
        ϕhguess= (ω_h_c.^σ).*Ycguess + ((ω_h_k.*Pkguess).^σ).*Ykguess

    else
        ϕhguess = [p_o_t.ϕht[1] exp.(Guess[1:N-2])' p_o_t.ϕht[N]]

    end

    #Measure of firms
    M0 = rt[1].measure
    sim=(measure = KLS3_measure_trans(p,rt,N,M0),)

    if s.transition_AD == 1 # See comment in line 202 of main.jl.
    sim=merge(sim,(w=zeros(eltype(p_o_t.w),N,1),ξ=zeros(eltype(p_o_t.w),N,1),Pk=zeros(eltype(p_o_t.w),N,1),ϕ_h=zeros(eltype(p_o_t.w),N,1),share_x=zeros(eltype(p_o_t.w),N,1),share_d=zeros(eltype(p_o_t.w),N,1),K=zeros(eltype(p_o_t.w),N,1),inv_agg=zeros(eltype(p_o_t.w),N,1),I=zeros(eltype(p_o_t.w),N,1),M=zeros(eltype(p_o_t.w),N,1),Yk=zeros(eltype(p_o_t.w),N,1),C=zeros(eltype(p_o_t.w),N,1),Yc=zeros(eltype(p_o_t.w),N,1),PdYd=zeros(eltype(p_o_t.w),N,1),PxYx=zeros(eltype(p_o_t.w),N,1),PxYx_USD=zeros(eltype(p_o_t.w),N,1),PmYm=zeros(eltype(p_o_t.w),N,1),tariffsincome=zeros(eltype(p_o_t.w),N,1),FC=zeros(eltype(p_o_t.w),N,1),N=zeros(eltype(p_o_t.w),N,1),n_supply=zeros(eltype(p_o_t.w),N,1),n_demand=zeros(eltype(p_o_t.w),N,1),mc_n=zeros(eltype(p_o_t.w),N,1),a_supply=zeros(eltype(p_o_t.w),N,1),a_demand=zeros(eltype(p_o_t.w),N,1),mc_a=zeros(eltype(p_o_t.w),N,1),y_supply=zeros(eltype(p_o_t.w),N,1),y_demand=zeros(eltype(p_o_t.w),N,1),mc_y=zeros(eltype(p_o_t.w),N,1),k_supply=zeros(eltype(p_o_t.w),N,1),k_demand=zeros(eltype(p_o_t.w),N,1),mc_k=zeros(eltype(p_o_t.w),N,1),mc_Yk_belief=zeros(eltype(p_o_t.w),N,1),mc_tariffs_belief=zeros(eltype(p_o_t.w),N,1),mc_y_belief=zeros(eltype(p_o_t.w),N,1),Sales=zeros(eltype(p_o_t.w),N,1),GDP=zeros(eltype(p_o_t.w),N,1),X_GDP=zeros(eltype(p_o_t.w),N,1),D_GDP=zeros(eltype(p_o_t.w),N,1),X_Sales=zeros(eltype(p_o_t.w),N,1),D_Sales=zeros(eltype(p_o_t.w),N,1),NX_GDP=zeros(eltype(p_o_t.w),N,1),NX_Sales=zeros(eltype(p_o_t.w),N,1),X_D=zeros(eltype(p_o_t.w),N,1),credit=zeros(eltype(p_o_t.w),N,1),credit_gdp=zeros(eltype(p_o_t.w),N,1),d_agg=zeros(eltype(p_o_t.w),N,1),NFA_GDP=zeros(eltype(p_o_t.w),N,1),k_wagebill=zeros(eltype(p_o_t.w),N,1),x_share_av=zeros(eltype(p_o_t.w),N,1),ln_sales_sd=zeros(eltype(p_o_t.w),N,1),ln_sales_d_sd=zeros(eltype(p_o_t.w),N,1),sales_sd=zeros(eltype(p_o_t.w),N,1),sales_d_sd=zeros(eltype(p_o_t.w),N,1),sales_avg_nx=zeros(eltype(p_o_t.w),N,1),labor_avg_nx=zeros(eltype(p_o_t.w),N,1),sales_d_avg_nx=zeros(eltype(p_o_t.w),N,1),sales_avg_x=zeros(eltype(p_o_t.w),N,1),labor_avg_x=zeros(eltype(p_o_t.w),N,1),labor_tot=zeros(eltype(p_o_t.w),N,1),labor_tot_x=zeros(eltype(p_o_t.w),N,1),labor_tot_nx=zeros(eltype(p_o_t.w),N,1),labor_tot_w=zeros(eltype(p_o_t.w),N,1),labor_x_share=zeros(eltype(p_o_t.w),N,1),labor_nx_share=zeros(eltype(p_o_t.w),N,1),sales_d_avg_x=zeros(eltype(p_o_t.w),N,1),xpremium_sales=zeros(eltype(p_o_t.w),N,1),xpremium_labor=zeros(eltype(p_o_t.w),N,1),xpremium_sales_d=zeros(eltype(p_o_t.w),N,1),ln_sales_sd_mean=zeros(eltype(p_o_t.w),N,1),ln_sales_d_sd_mean=zeros(eltype(p_o_t.w),N,1),sales_sd_mean=zeros(eltype(p_o_t.w),N,1),sales_d_sd_mean=zeros(eltype(p_o_t.w),N,1),labor_sd_mean=zeros(eltype(p_o_t.w),N,1),avg_productivity1=zeros(eltype(p_o_t.w),N,1),avg_productivity2=zeros(eltype(p_o_t.w),N,1),avg_productivity3=zeros(eltype(p_o_t.w),N,1),X_Laspeyres=zeros(eltype(p_o_t.w),N,1),D_Laspeyres=zeros(eltype(p_o_t.w),N,1),Sales_Laspeyres=zeros(eltype(p_o_t.w),N,1),GDP_Laspeyres=zeros(eltype(p_o_t.w),N,1),Imports=zeros(eltype(p_o_t.w),N,1),Imports_C=zeros(eltype(p_o_t.w),N,1),Imports_K=zeros(eltype(p_o_t.w),N,1),Imports_C_Laspeyres=zeros(eltype(p_o_t.w),N,1),Imports_K_Laspeyres=zeros(eltype(p_o_t.w),N,1),Imports_Laspeyres=zeros(eltype(p_o_t.w),N,1),a_min_share=zeros(eltype(p_o_t.w),N,1),a_max_share=zeros(eltype(p_o_t.w),N,1)))
    else
    sim=merge(sim,(w=zeros(N,1),ξ=zeros(N,1),Pk=zeros(N,1),ϕ_h=zeros(N,1),share_x=zeros(N,1),share_d=zeros(N,1),K=zeros(N,1),inv_agg=zeros(N,1),I=zeros(N,1),M=zeros(N,1),Yk=zeros(N,1),C=zeros(N,1),Yc=zeros(N,1),PdYd=zeros(N,1),PxYx=zeros(N,1),PxYx_USD=zeros(N,1),PmYm=zeros(N,1),tariffsincome=zeros(N,1),FC=zeros(N,1),N=zeros(N,1),n_supply=zeros(N,1),n_demand=zeros(N,1),mc_n=zeros(N,1),a_supply=zeros(N,1),a_demand=zeros(N,1),mc_a=zeros(N,1),y_supply=zeros(N,1),y_demand=zeros(N,1),mc_y=zeros(N,1),k_supply=zeros(N,1),k_demand=zeros(N,1),mc_k=zeros(N,1),mc_Yk_belief=zeros(N,1),mc_tariffs_belief=zeros(N,1),mc_y_belief=zeros(N,1),Sales=zeros(N,1),GDP=zeros(N,1),X_GDP=zeros(N,1),D_GDP=zeros(N,1),X_Sales=zeros(N,1),D_Sales=zeros(N,1),NX_GDP=zeros(N,1),NX_Sales=zeros(N,1),X_D=zeros(N,1),credit=zeros(N,1),credit_gdp=zeros(N,1),d_agg=zeros(N,1),NFA_GDP=zeros(N,1),k_wagebill=zeros(N,1),x_share_av=zeros(N,1),ln_sales_sd=zeros(N,1),ln_sales_d_sd=zeros(N,1),sales_sd=zeros(N,1),sales_d_sd=zeros(N,1),sales_avg_nx=zeros(N,1),labor_avg_nx=zeros(N,1),sales_d_avg_nx=zeros(N,1),sales_avg_x=zeros(N,1),labor_avg_x=zeros(N,1),labor_tot=zeros(N,1),labor_tot_x=zeros(N,1),labor_tot_nx=zeros(N,1),labor_tot_w=zeros(N,1),labor_x_share=zeros(N,1),labor_nx_share=zeros(N,1),sales_d_avg_x=zeros(N,1),xpremium_sales=zeros(N,1),xpremium_labor=zeros(N,1),xpremium_sales_d=zeros(N,1),ln_sales_sd_mean=zeros(N,1),ln_sales_d_sd_mean=zeros(N,1),sales_sd_mean=zeros(N,1),sales_d_sd_mean=zeros(N,1),labor_sd_mean=zeros(N,1),avg_productivity1=zeros(N,1),avg_productivity2=zeros(N,1),avg_productivity3=zeros(N,1),X_Laspeyres=zeros(N,1),D_Laspeyres=zeros(N,1),Sales_Laspeyres=zeros(N,1),GDP_Laspeyres=zeros(N,1),Imports=zeros(N,1),Imports_C=zeros(N,1),Imports_K=zeros(N,1),Imports_C_Laspeyres=zeros(N,1),Imports_K_Laspeyres=zeros(N,1),Imports_Laspeyres=zeros(N,1),a_min_share=zeros(N,1),a_max_share=zeros(N,1)))
    end

    for t=2:N
    ### Shocks

        p_o_t=merge(p_o_t,(r = rv[t], # shock to interest rate
        r_lag = rv[t-1],
        Pf_lag = Pfv[t-1], # shock to foreign price
        Pf = Pfv[t],
        Yf = Yfv[t], # shock to foreign demand
        θ = θ_v[t], # shock to collateral constraint
        β = β_v[t],    # shock to discount factor
        δ = δ_v[t], # shock to depreciation rate
        τ = τ_v[t], # shock to iceberg costs
        τ_x = τ_x_v[t], # shocks to tariffs
        τ_m_c = τ_m_c_v[t], # shocks to consumption tariffs
        τ_m_k = τ_m_k_v[t], # shocks to capital tariffs

        ## GE prices
        ϕ_h = ϕhguess[t],
        w = wguess[t],
        ξ = ξguess[t],
        Pk = Pkguess[t],
        Pk_lag = Pkguess[t-1]))

        if s.tariffsincome==1
            p_o_t=merge(p_o_t,(Yk = Ykguess[t],
            Yc = Ycguess[t],
            tariffsincome=p_o_t.tariffsincomet[t]))
        end

        ## Useful objects

        n = rt[t].n
        mat = rt[t].m  # intermediates
        measure = sim.measure[:,:,t]  # measure
        pd = rt[t].pd
        yd = rt[t].yd
        pf = rt[t].pf
        yf = rt[t].yf
        e = rt[t].e
        k = rt[t].k
        if t<N
            kp = rt[t+1].k
        end
        F_mat = rt[t].F_mat
        c = rt[t].c

        # Saving prices

        sim.w[t,1]=p_o_t.w
        sim.ξ[t,1]=p_o_t.ξ
        sim.Pk[t,1] = p_o_t.Pk
        sim.ϕ_h[t,1] = p_o_t.ϕ_h

        # Exporters and non-exporters
        sim.share_x[t,1] = sum(measure.*e) # share of exporters
        sim.share_d[t,1] = sum(measure.*(1 .-e))  #share of non-exporters

        ## Capital good
        sim.K[t,1] = sum(real(measure.*k))
        if t>1 && t<N
            measurep =  sim.measure[:,:,t+1]
            sim.inv_agg[t,1] = real(sum(measurep.*kp) - (1 .-δ_v[t])*sum(measure.*k))
        elseif t==N #For last period (t=N)
            # We assume that we reached steady state when t=N
            sim.inv_agg[t,1] = δ*sim.K[t,1]
        end

        sim.I[t,1] = max(sim.inv_agg[t,1],0) #Demand for new investment goods
        sim.M[t,1] = sum(real(measure.*mat))

        yd_k = ((pd/(p_o_t.Pk*ω_h_k)).^(-σ)) * (sim.I[t,1] + sim.M[t,1])
        ym_k = ((p_o_t.ξ*Pm_k*(1+τ_m_k)/(p_o_t.Pk*ω_m_k)).^(-σ)) * (sim.I[t,1] + sim.M[t,1])

        sim.Yk[t,1] = (sum(real(measure.*ω_h_k.*(yd_k.^((σ-1)/σ)))) + ω_m_k*(ym_k.^((σ-1)/σ)))^(σ/(σ-1))
        sim.Pk[t,1] = (sum(real(measure.*(ω_h_k^σ).*(pd.^(1-σ)))) + (ω_m_k^σ)*((p_o_t.ξ*Pm_k*(1+τ_m_k)).^(1-σ)) )^(1/(1-σ))

        ## Consumption good

        sim.C[t,1] = sum(real(measure.*c))

        #yd_c = ((pd/ω_h_c).^(-σ)) * sim.C[t,1]
        yd_c = max.(real(yd - yd_k),0.00001)
        ym_c = (p_o_t.ξ*Pm_c*(1+τ_m_c)/(ω_m_c)).^(-σ) * sim.C[t,1]
        sim.Yc[t,1] = (sum(real(measure.*ω_h_c.*(yd_c.^((σ-1)/σ)))) + ω_m_c*(ym_c.^((σ-1)/σ)) )^(σ/(σ-1))

        ## ϕ_h

        sim.ϕ_h[t,1] = (ω_h_c^σ)*sim.Yc[t,1] + ((ω_h_k*p_o_t.Pk)^σ)*sim.Yk[t,1]

        ## Compute price and quantity indexes

        #Domestic sales
         sim.PdYd[t,1] =sum(real(measure[real(pd).>0].*(pd[real(pd).>0].*yd[real(pd).>0])))

        #Exports
         sim.PxYx[t,1] = p_o_t.ξ*sum(real(measure[real(pf).>0].*(pf[real(pf).>0].*yf[real(pf).>0])))

         sim.PxYx_USD[t,1] = sim.PxYx[t,1]/p_o_t.ξ

         #Imports
        sim.PmYm[t,1] = p_o_t.ξ*(1+τ_m_c)*Pm_c*ym_c + p_o_t.ξ*(1+τ_m_k)*Pm_k*ym_k

        # Tariffs Income (T)
        sim.tariffsincome[t,1]=p_o_t.ξ*τ_m_c*Pm_c*ym_c + p_o_t.ξ*τ_m_k*Pm_k*ym_k
        if s.tariffsincome_PE==1
            sim.tariffsincome[t,1] = p_o_t.tariffsincome
        end
        ## Compute aggregate variables needed to evaluate market clearing conditions

        #Labor and capital
        sim.FC[t,1] =  sum(measure.*e.*F_mat)
        sim.N[t,1] = sum(real(measure.*n))    +((sim.FC[t,1])/p_o_t.w)*(1 .-s.fcost_fgoods)

        ## Market clearing conditions

        #Labor
        sim.n_supply[t,1] = 1
        sim.n_demand[t,1] = sim.N[t,1]
        #sim.mc_n[t,1] = ((1+sim.n_demand[t,1])/(1+sim.n_supply[t,1]))-1
        sim.mc_n[t,1] = log.(sim.n_demand[t,1]/sim.n_supply[t,1])

        #Assets
        #This market clearing condition is the result of reformulating the debt market clearing condition
        #In the original model, the sum of debt has to equal zero. Once the model is reformulated, this condition becomes the one below.
        sim.a_supply[t,1] = sum(real(measure.*rt[t].a))
        sim.a_demand[t,1] = sim.K[t,1]
        #sim.mc_a[t,1] = (1+sim.a_demand[t,1])/(1+sim.a_supply[t,1])-1

        sim.mc_a[t,1] = log.(sim.a_demand[t,1]/sim.a_supply[t,1])

        #Final goods
        sim.y_supply[t,1] = sim.Yc[t,1]
        sim.y_demand[t,1] = sim.C[t,1]
        sim.mc_y[t,1] = log.(sim.y_demand[t,1]/sim.y_supply[t,1])

        # Capital goods
        sim.k_supply[t,1] = sim.Yk[t,1]
        sim.k_demand[t,1] = sim.I[t,1] + sim.M[t,1] + sim.FC[t,1]*s.fcost_fgoods
        sim.mc_k[t,1] = log.(sim.k_demand[t,1]/sim.k_supply[t,1])
#       sim.mc_k[t,1] =sim.k_demand[t,1]/sim.k_supply[t,1]-1

        #Beliefs
        #To solve the entrepreneur's problem, they need to have a belief about the aggregate price and quantity indexes in the market
        #In equilibrium, these beliefs need to be consistent with the actual aggregate prices and quantities
        #sim.mc_y_belief[t,1] = log.(sim.Yh[t,1]/m.Yh) sim.mc_y_belief = 10*log.((1+sim.Ycpi_t)/(1+m.Yh))  sim.mc_y_belief[t,1] = ((1+sim.Yh[t,1])/(1+m.Yh))-1

        if s.tariffsincome==1
            sim.mc_Yk_belief[t,1] = log.(sim.Yk[t,1]/p_o_t.Yk)
            sim.mc_tariffs_belief[t,1] = log.(sim.tariffsincome[t,1]/p_o_t.tariffsincome)
            sim.mc_y_belief[t,1] = ((1+sim.Yc[t,1])/(1+p_o_t.Yc))-1
        else
            sim.mc_Yk_belief[t,1] = 0
            sim.mc_tariffs_belief[t,1] = 0
            sim.mc_y_belief[t,1] = log.(sim.ϕ_h[t,1]/ϕhguess[t])
        end

        ## Main statistics

        sim.Sales[t,1] = sim.PdYd[t,1] + sim.PxYx[t,1] #sim.PxYx is already denominated in units of the domestic final good
        sim.GDP[t,1] = sim.Sales[t,1] -  sim.M[t,1]*p_o_t.Pk #

        sim.X_GDP[t,1] = sim.PxYx[t,1]/sim.GDP[t,1]
        sim.D_GDP[t,1]  = sim.PdYd[t,1] /sim.GDP[t,1]

        sim.X_Sales[t,1] = sim.PxYx[t,1]/sim.Sales[t,1]
        sim.D_Sales[t,1] = sim.PdYd[t,1]/sim.Sales[t,1]

         sim.NX_GDP[t,1] = (sim.PxYx[t,1]-sim.PmYm[t,1])/sim.GDP[t,1]
         sim.NX_Sales[t,1] = (sim.PxYx[t,1]-sim.PmYm[t,1])/sim.Sales[t,1]

         sim.X_D[t,1] = sim.PxYx[t,1]/sim.PdYd[t,1]

        # Every credit statistic only for tradables
        rt[t]=merge(rt[t],(d = (1+r).*(p_o_t.Pk_lag*k - rt[t].a),))

        sim.credit[t,1] = sum(real(measure.*max.(real(rt[t].d),0))) # only entrepreneurs/firms (for workers r.d is negative)
        sim.credit_gdp[t,1] = sim.credit[t,1]/sim.GDP[t,1]
        sim.d_agg[t,1] = sum(real(measure.*rt[t].d)) # for both firms and workers
        sim.NFA_GDP[t,1] = -sim.d_agg[t,1]/sim.GDP[t,1]

        sim.k_wagebill[t,1] = p_o_t.Pk_lag*sim.K[t,1]/(p_o_t.w*sim.n_supply[t,1])

        # Sales
        sales_x =  p_o_t.ξ*rt[t].pf.*rt[t].yf #r.pf is denominated in foreign currency, so we adjust it
        sales_d = rt[t].pd.*rt[t].yd
        sales = sales_d+sales_x
        x_share = sales_x./sales

        sim.x_share_av[t,1] =  sum(filter(!isnan,real(measure.*e.* x_share))) / sum(real(measure.*e))
        # sim.x_share_av[isnan.(sim.x_share_av[t,1])].=1 #Don't know the role of this, but it throws error in Julia given sim.x_share_av is a scalar

        if s.extra_results==1
            ln_sales = log.(sales)
            ln_sales_d =log.(sales_d)
            ln_sales_ind = isnan.(ln_sales).==0
            ln_sales_d_ind = isnan.(ln_sales_d).==0

            ln_sales_mean = sum(real(measure[ln_sales_ind].*ln_sales[ln_sales_ind]))
            sim.ln_sales_sd[t,1] = sqrt.(sum(real(measure[ln_sales_ind] .* (ln_sales[ln_sales_ind] .- ln_sales_mean).^2)))

            ln_sales_d_mean = sum(real(measure[ln_sales_d_ind] .* ln_sales_d[ln_sales_d_ind]))
            sim.ln_sales_d_sd[t,1] = sqrt.( sum(real(measure[ln_sales_d_ind] .* (ln_sales_d[ln_sales_d_ind] .- ln_sales_d_mean).^2)))

            sales_mean = sum(real(measure.*sales))
            sim.sales_sd[t,1] = sqrt.(sum(real(measure .* (sales .- sales_mean).^2)))

            sales_d_mean = sum(real(measure .* sales_d))
            sim.sales_d_sd[t,1] = sqrt.(sum(real(measure .* (sales_d .- sales_d_mean).^2)))

            sim.sales_avg_nx[t,1] = sum(real(measure.*(1 .-e).*sales)) / sum(real(measure.*(1 .-e)))
            sim.labor_avg_nx[t,1] = sum(real(measure.*(1 .-e).*(n))) / sum(real(measure.*(1 .-e)))
            sim.sales_d_avg_nx[t,1] = sim.sales_avg_nx[t,1]

            sim.sales_avg_x[t,1] = sum(real(measure.*e.*sales))/ sum(real(measure.*e))
            sim.labor_avg_x[t,1] = ((sim.FC[t,1]/p_o_t.w)*(1-s.fcost_fgoods) + sum(real(measure.*e.*(n)))) / sum(real(measure.*e))

        # These include fixed costs
            sim.labor_tot[t,1] = sim.N[t,1]  # total labor demand for production
            sim.labor_tot_x[t,1] = (sim.FC[t,1]/p_o_t.w)*(1-s.fcost_fgoods)  + sum(real(measure.*e.*(n)))
            sim.labor_tot_nx[t,1] =  sum(real(measure.*(1 .-e).*(n)))
            sim.labor_tot_w[t,1] = sum(real(measure))

            sim.labor_x_share[t,1] = sim.labor_tot_x[t,1]/sim.labor_tot[t,1]
            sim.labor_nx_share[t,1] = sim.labor_tot_nx[t,1]/sim.labor_tot[t,1]

            sim.sales_d_avg_x[t,1] = sum(real(measure.*e.*sales_d))/ sum(real(measure.*e))

            sim.xpremium_sales[t,1] = sim.sales_avg_x[t,1]/sim.sales_avg_nx[t,1]
            sim.xpremium_labor[t,1] = sim.labor_avg_x[t,1]/sim.labor_avg_nx[t,1]

            sim.xpremium_sales_d[t,1] = sim.sales_d_avg_x[t,1]/sim.sales_d_avg_nx[t,1]

            sim.ln_sales_sd_mean[t,1] = sim.ln_sales_sd[t,1] /  ln_sales_mean
            sim.ln_sales_d_sd_mean[t,1] = sim.ln_sales_d_sd[t,1] /  ln_sales_d_mean
            sim.sales_sd_mean[t,1] = sim.sales_sd[t,1] /  sales_mean
            sim.sales_d_sd_mean[t,1] = sim.sales_d_sd[t,1] /  sales_d_mean

            labor_mean= sum(real(measure.*(n+ (e.*F_mat./p_o_t.w )*(1-s.fcost_fgoods))))
            labor_sd=sqrt.(sum(real(measure.*( n+ (e.*F_mat./p_o_t.w)*(1-s.fcost_fgoods) .- labor_mean ).^2)))
            sim.labor_sd_mean[t,1] = labor_sd / labor_mean

           # Average productivity
            sim.avg_productivity1[t,1] = ((sum(real(measure.*(rt[t].z).^(σ-1))))/sum(measure) ).^(1/(σ - 1))
            sim.avg_productivity2[t,1] =  sum(real(measure.*sales.* rt[t].z)) / sum(real(measure .* sales))
            sim.avg_productivity3[t,1] =  sum(real(measure.*sales_d.* rt[t].z)) / sum(real(measure .* sales_d ))
        end

    ## Real Statistics
    #Real GDP, Real Exports and Real Domestic Sales

        sim.X_Laspeyres[t,1] = sum(real(measure.*p_o_t.ξt[1].*rt[1].pf.*rt[t].yf))
        sim.D_Laspeyres[t,1] = sum(real(measure.*rt[1].pd.*rt[t].yd))
        sim.Sales_Laspeyres[t,1] = sum(real(measure.*(rt[1].pd.*rt[t].yd+p_o_t.ξt[1].*rt[1].pf.*rt[t].yf)))
        sim.GDP_Laspeyres[t,1] = sum(real(measure.*(rt[1].pd.*rt[t].yd+p_o_t.ξt[1]*rt[1].pf.*rt[t].yf-mat*p_o_t.Pkt[1])))
        sim.Imports[t,1] = sim.PmYm[t,1]
        sim.Imports_C[t,1] = p_o_t.ξ*(1+τ_m_c)*Pm_c*ym_c
        sim.Imports_K[t,1] =  p_o_t.ξ*(1+τ_m_k)*Pm_k*ym_k
        sim.Imports_C_Laspeyres[t,1] = p_o_t.ξt[1].*(1+τ_m_c_v[1])*Pm_c*ym_c
        sim.Imports_K_Laspeyres[t,1] =  p_o_t.ξt[1].*(1+τ_m_k_v[1])*Pm_k*ym_k
        sim.Imports_Laspeyres[t,1] =sim.Imports_C_Laspeyres[t,1] +sim.Imports_K_Laspeyres[t,1]

     ## Solution-related statistics

        sim.a_min_share[t,1] = sum(real(measure[1,:]))
        sim.a_max_share[t,1] = sum(real(measure[end,:]))

    end
return sim,rt
end

############################# Transition_Vec2 (2 versions, first for using NLsolve [in-place function], second standard)#############################

####### For NLsolve #######
#This version of KLS3_transition_vec2 is in-place, that is, its output (market clearing conditions) is a modified value of an input (F, in this case). This is done because NLsolve needs this format to work. It is customary to add a ! at the end of in-place functions.

function KLS3_transition_vec2!(F,Guess,p,p_o_t,rt)
@unpack N, ω_h_c, σ, ω_h_k,rv,Pfv,Yfv,θ_v,β_v,δ_v,τ_v,τ_x_v,τ_m_c_v,τ_m_k_v,F_base,ω_m_k,ω_m_c,Pm_c,τ_m_c,τ_m_k,tariffsincome,Pm_k = p #See comment at the beginning of KL3_measure_trans

##Transform vector of guess into standard matrix form (used for LeastSquaresOptim solver)
if s.solver_LM==1
Guess1=copy(Guess)
Guess=hcat(Guess1[1:N-2],Guess1[N-1:2*N-4],Guess1[2*N-3:3*N-6],Guess1[3*N-5:4*N-8],Guess1[4*N-7:5*N-10])
end

	wguess = [p_o_t.wt[1] exp.(Guess[N-1:2N-4])' p_o_t.wt[N]]
	ξguess = [p_o_t.ξt[1] exp.(Guess[2N-3:3N-6])' p_o_t.ξt[N]]
	Pkguess = [p_o_t.Pkt[1] exp.(Guess[3N-5:4N-8])' p_o_t.Pkt[N]]

    if s.tariffsincome == 1
		Ycguess = [p_o_t.Yct[1] exp.(Guess[1:N-2])' p_o_t.Yct[N]]
        Ykguess = [p_o_t.Ykt[1] exp.(Guess[4*N-7:5*N-10])' p_o_t.Ykt[N]]

        ϕhguess= (ω_h_c.^σ).*Ycguess + ((ω_h_k.*Pkguess).^σ).*Ykguess

    else
		ϕhguess = [p_o_t.ϕht[1] exp.(Guess[1:N-2])' p_o_t.ϕht[N]]

    end

    for t=N-1:-1:2

        ## Shocks

		# shock to interest rate
        p_o_t=merge(p_o_t,(r = p.rv[t],
        r_lag = p.rv[t-1],

        # shock to foreign price
        Pf_lag = p.Pfv[t-1],
        Pf = p.Pfv[t],

        # shock to foreign demand
        Yf = p.Yfv[t],

        # shock to collateral constraint
        θ = p.θ_v[t],

        # shock to discount factor
        β = p.β_v[t],

        # shock to depreciation rate
        δ = p.δ_v[t],

        # shock to iceberg costs
        τ = p.τ_v[t],

        # shock to tariffs
        τ_x = p.τ_x_v[t],
        τ_m_c = p.τ_m_c_v[t],
        τ_m_k = p.τ_m_k_v[t],

        ## GE prices ##

        ϕ_h = ϕhguess[t],
        ϕ_h_lag = ϕhguess[t-1],

        w = wguess[t],
        w_lag = wguess[t-1],

        ξ = ξguess[t],
        ξ_lag = ξguess[t-1],

        Pk = Pkguess[t],
        Pk_lag = Pkguess[t-1]))

		## Fixed costs
        if  s.fcost_fgoods==0 # If in units of labor
            p_o_t=merge(p_o_t,(F = p_o_t.w*F_base,))
         else
            p_o_t=merge(p_o_t,(F = F_base,))
        end

      ## Tariff income (initialized to 0) ##
        if s.tariffsincome == 1

            p_o_t=merge(p_o_t,(Yk = Ykguess[t],
            Yc = Ycguess[t]))
            #Yc=p_o_t.Yc
            #Yc =  (ϕhguess[t] - (ω_m_k*Pkguess[t])^σ*Ykguess[t])./(ω_m_c^σ)
            ym_c = p_o_t.Yc*(p_o_t.ξ*Pm_c*(1+τ_m_c)/ω_m_c)^(-σ)
            ym_k = p_o_t.Yk*(Pkguess[t]^σ)*(p_o_t.ξ*Pm_k*(1+τ_m_k)/ω_m_k)^(-σ)
            if s.tariffsincome_PE==0
                p_o_t=merge(p_o_t,(tariffsincome = τ_m_c*p_o_t.ξ*Pm_c*ym_c + τ_m_k*p_o_t.ξ*Pm_k*ym_k,))
			else
				p_o_t=merge(p_o_t,(tariffsincome = tariffsincome,))
            end
            p_o_t.tariffsincomet[t]=p_o_t.tariffsincome
        end

        ## Solve static and dynamic problems
        # Period 2
		if t==2
			p_o_temp = KLS3_staticproblem_period2_altTiming(p_o_t,p)
		else
			p_o_temp = KLS3_staticproblem(p_o_t,p)
		end

        #Fixed costs
		p_o_temp=merge(p_o_temp,(S_mat = p_o_temp.e*p_o_t.F,
        F_mat = p_o_temp.e*p_o_t.F,
        profits=p_o_t.w .+ p_o_t.tariffsincome .+ p_o_temp.e.*p_o_temp.π_x .+ (1 .-p_o_temp.e).*p_o_temp.π_nx))
        rt[t]=p_o_temp

    end

	rt = KLS3_dynamicproblem_trans_vec_t(p_o_t,p,rt)

	if s.flag_simulate == 0
		sim_fun, rt = KLS3_simulate_trans(p,p_o_t,rt,Guess)
	end

    # Market clearing conditions

    if s.tariffsincome==0
        for j=1:N-2
        F[j]=sim_fun.mc_n[j+1]
        end
        for j=N-1:2*N-4
        F[j]=sim_fun.mc_y[j+1-N+2]
        end
        for j=2*N-3:3*N-6
        F[j]=sim_fun.mc_k[2:N-1]
        end
        for j=3*N-5:4*N-8
        F[j]=sim_fun.mc_y_belief[2:N-1]
        end
    elseif s.tariffsincome==1
        for j=1:N-2
        F[j]=sim_fun.mc_n[j+1]
        end
        for j=N-1:2*N-4
        F[j]=sim_fun.mc_y[j+3-N]
        end
        for j=2*N-3:3*N-6
        F[j]=sim_fun.mc_k[j+5-2*N]
        end
        for j=3*N-5:4*N-8
        F[j]=sim_fun.mc_y_belief[j+7-3*N]
        end
        for j=4*N-7:5*N-10
        F[j]=sim_fun.mc_Yk_belief[j+9-4*N]
        end
    end
    return F
end

####################### Standard #######################

function KLS3_transition_vec2(Guess,p,p_o_t,rt)
#See comments of the in-place version (KLS3_transition_vec2!)

@unpack N, ω_h_c, σ, ω_h_k,rv,Pfv,Yfv,θ_v,β_v,δ_v,τ_v,τ_x_v,τ_m_c_v,τ_m_k_v,F_base,ω_m_k,ω_m_c,Pm_c,τ_m_c,τ_m_k,tariffsincome,Pm_k=p

    wguess = [p_o_t.wt[1] exp.(Guess[N-1:2N-4])' p_o_t.wt[N]]
    ξguess = [p_o_t.ξt[1] exp.(Guess[2N-3:3N-6])' p_o_t.ξt[N]]
    Pkguess = [p_o_t.Pkt[1] exp.(Guess[3N-5:4N-8])' p_o_t.Pkt[N]]

    if s.tariffsincome == 1
        Ycguess = [p_o_t.Yct[1] exp.(Guess[1:N-2])' p_o_t.Yct[N]]
        Ykguess = [p_o_t.Ykt[1] exp.(Guess[4*N-7:5*N-10])' p_o_t.Ykt[N]]

        ϕhguess= (ω_h_c.^σ).*Ycguess + ((ω_h_k.*Pkguess).^σ).*Ykguess

    else
        ϕhguess = [p_o_t.ϕht[1] exp.(Guess[1:N-2])' p_o_t.ϕht[N]]

    end

    for t=N-1:-1:2

        ## Shocks

        #global p_o_t

        # shock to interest rate
        p_o_t=merge(p_o_t,(r = p.rv[t],
        r_lag = p.rv[t-1],

        # shock to foreign price
        Pf_lag = p.Pfv[t-1],
        Pf = p.Pfv[t],

        # shock to foreign demand
        Yf = p.Yfv[t],

        # shock to collateral constraint
        θ = p.θ_v[t],

        # shock to discount factor
        β = p.β_v[t],

        # shock to depreciation rate
        δ = p.δ_v[t],

        # shock to iceberg costs
        τ = p.τ_v[t],

        # shock to tariffs
        τ_x = p.τ_x_v[t],
        τ_m_c = p.τ_m_c_v[t],
        τ_m_k = p.τ_m_k_v[t],

        ## GE prices ##

        ϕ_h = ϕhguess[t],
        ϕ_h_lag = ϕhguess[t-1],

        w = wguess[t],
        w_lag = wguess[t-1],

        ξ = ξguess[t],
        ξ_lag = ξguess[t-1],

        Pk = Pkguess[t],
        Pk_lag = Pkguess[t-1]))

        ## Fixed costs $$
        if  s.fcost_fgoods==0 # If in units of labor
            p_o_t=merge(p_o_t,(F = p_o_t.w*F_base,))
         else
            p_o_t=merge(p_o_t,(F = F_base,))
        end

      ## Tariff income (initialized to 0) ##
        if s.tariffsincome == 1

            p_o_t=merge(p_o_t,(Yk = Ykguess[t],
            Yc = Ycguess[t]))
            #Yc=p_o_t.Yc
            #Yc =  (ϕhguess[t] - (ω_m_k*Pkguess[t])^σ*Ykguess[t])./(ω_m_c^σ)
            ym_c = p_o_t.Yc*(p_o_t.ξ*Pm_c*(1+τ_m_c)/ω_m_c)^(-σ)
            ym_k = p_o_t.Yk*(Pkguess[t]^σ)*(p_o_t.ξ*Pm_k*(1+τ_m_k)/ω_m_k)^(-σ)
            if s.tariffsincome_PE==0
                p_o_t=merge(p_o_t,(tariffsincome = τ_m_c*p_o_t.ξ*Pm_c*ym_c + τ_m_k*p_o_t.ξ*Pm_k*ym_k,))
			else
				p_o_t=merge(p_o_t,(tariffsincome = tariffsincome,))
            end
            p_o_t.tariffsincomet[t]=p_o_t.tariffsincome
        end

        ## Solve static and dynamic problems ##

        # Period 2
        if t==2
            p_o_temp = KLS3_staticproblem_period2_altTiming(p_o_t,p)
        else
            p_o_temp = KLS3_staticproblem(p_o_t,p)
        end

        #Fixed costs
        p_o_temp=merge(p_o_temp,(S_mat = p_o_temp.e*p_o_t.F,
        F_mat = p_o_temp.e*p_o_t.F,
        profits=p_o_t.w .+ p_o_t.tariffsincome .+ p_o_temp.e.*p_o_temp.π_x .+ (1 .-p_o_temp.e).*p_o_temp.π_nx))
        rt[t]=p_o_temp
    end

    rt = KLS3_dynamicproblem_trans_vec_t(p_o_t,p,rt)

    #ap_ind_tr_Julia=zeros(100,100,20)
    #for t=1:20
        #ap_ind_tr_Julia[:,:,t]=rt[t].ap_ind
    #end
    #file=matopen("ap_ind_tr_Julia.mat","w")
    #write(file,"ap_ind_tr_Julia",ap_ind_tr_Julia)
    #close(file)

    if s.flag_simulate == 0
        sim_fun, rt = KLS3_simulate_trans(p,p_o_t,rt,Guess)
    end

    # Market clearing conditions

    if s.tariffsincome==0
         mc = [sim_fun.mc_n[2:N-1]' sim_fun.mc_y[2:N-1]' sim_fun.mc_k[2:N-1]' sim_fun.mc_y_belief[2:N-1]']'
    elseif s.tariffsincome==1
        mc = [sim_fun.mc_n[2:N-1]' sim_fun.mc_y[2:N-1]' sim_fun.mc_k[2:N-1]' sim_fun.mc_y_belief[2:N-1]' sim_fun.mc_Yk_belief[2:N-1]']'
    end

    return mc, p_o_t, sim_fun, rt
end


text:
import Compat: @static, is_windows, is_apple

@static if is_windows()
    error("The shoco C library doesn't support Windows")
end

using BinDeps

@BinDeps.setup

shoco = library_dependency("shoco")

if isdir(srcdir(shoco))
    rm(srcdir(shoco), recursive=true)
    mkdir(srcdir(shoco))
end

if isdir(BinDeps.downloadsdir(shoco))
    rm(BinDeps.downloadsdir(shoco), recursive=true)
    mkdir(BinDeps.downloadsdir(shoco))
end

sha = "4dee0fc850cdec2bdb911093fe0a6a56e3623b71"

provides(Sources, URI("https://github.com/Ed-von-Schleck/shoco/archive/$(sha).zip"), shoco,
         unpacked_dir="shoco-$sha")

provides(BuildProcess, (@build_steps begin
    GetSources(shoco)
    CreateDirectory(joinpath(BinDeps.builddir(shoco), "shoco"))
    @build_steps begin
        ChangeDirectory(joinpath(BinDeps.builddir(shoco), "shoco"))
        FileRule(joinpath(libdir(shoco), "shoco." * BinDeps.shlib_ext), @build_steps begin
            CreateDirectory(libdir(shoco))
            CCompile(joinpath(srcdir(shoco), "shoco-$sha", "shoco.c"),
                     joinpath(libdir(shoco), "shoco." * BinDeps.shlib_ext),
                     ["-fPIC", "-std=c99", is_apple() ? "-dynamiclib" : "-shared"], [])
        end)
    end
end), shoco)

@BinDeps.install Dict(:shoco => :shoco)

ext: jl


text:
module TidalFluxCalibrations

using Reexport

@reexport using TidalFluxQuantities, StatsBase, Interpolations

export CalibrationModel,
    PolynomialCalibrationModel,
    Calibration,
    to_quantity,
    from_quantity

abstract type CalibrationModel{T<:Quantity, F<:Quantity} <: RegressionModel end

include("calibrations.jl")

"""
    interpolateCalibration(c::Calibration)

Takes a Calibration and returns the from quantity
that has been interpolated to the times of the
to quantity
"""
function interpolatecal{T,F}(cal::Calibration{T,F})
    qi = interpolate((times(from_quantity(cal)),),
                     quantity(from_quantity(cal)),
                     Gridded(Linear())
                     )
    ts = times(to_quantity(cal))
    qs = qi[ts]
    T(ts,qs)
end

include("polynomial.jl")

end # module


text:
abstract type AbstractMultiThreading <: ParallelMode end
struct MultiThreading{abstractMode <: BaseMode} <: AbstractMultiThreading
    numberOfThreads::Int64
    sub_mod::abstractMode
    MultiThreading(sub_mod::abstractMode = SerialMode()) where {abstractMode <: BaseMode} = new{abstractMode}(Int64(Threads.nthreads()), sub_mod)
end

function pricer_macro_multithreading(model_type, payoff_type)
    @eval begin
        function pricer(mcProcess::$model_type, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration{<:Integer, <:Integer, <:AbstractMonteCarloMethod, <:AbstractMultiThreading, <:Random.AbstractRNG}, abstractPayoff::$payoff_type)
            zero_typed = predict_output_type_zero_(mcProcess, rfCurve, mcConfig, abstractPayoff)
            price = zeros(typeof(zero_typed), Threads.nthreads())
            Threads.@threads for i = 1:Threads.nthreads()
                price[i] = pricer(mcProcess, rfCurve, MonteCarloConfiguration(div(mcConfig.Nsim, Threads.nthreads()), mcConfig.Nstep, mcConfig.monteCarloMethod, mcConfig.parallelMode.sub_mod, mcConfig.seed), abstractPayoff)
            end
            Out = sum(price) / Threads.nthreads()
            return Out
        end
    end
end

pricer_macro_multithreading(BaseProcess, AbstractPayoff)
pricer_macro_multithreading(BaseProcess, Dict{AbstractPayoff, Number})
pricer_macro_multithreading(Dict{String, AbstractMonteCarloProcess}, Dict{String, Dict{AbstractPayoff, Number}})
pricer_macro_multithreading(VectorialMonteCarloProcess, Array{Dict{AbstractPayoff, Number}})

function pricer(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration{<:Integer, <:Integer, <:AbstractMonteCarloMethod, <:AbstractMultiThreading, <:Random.AbstractRNG}, abstractPayoffs::Array{abstractPayoff_}) where {abstractPayoff_ <: AbstractPayoff}
    zero_typed = predict_output_type_zero_(mcProcess, rfCurve, mcConfig, abstractPayoffs)
    price = zeros(typeof(zero_typed), length(abstractPayoffs), Threads.nthreads())
    Threads.@threads for i = 1:Threads.nthreads()
        price[:, i] = pricer(mcProcess, rfCurve, MonteCarloConfiguration(div(mcConfig.Nsim, Threads.nthreads()), mcConfig.Nstep, mcConfig.monteCarloMethod, mcConfig.parallelMode.sub_mod, mcConfig.seed), abstractPayoffs)
    end
    Out = sum(price, dims = 2) / Threads.nthreads()
    return Out
end


text:
using CSV, DataFrames, Distributed, Dates, LinearAlgebra, Distributions, DelimitedFiles, SharedArrays

# Custom package
using Jevo

# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

# Get number of workers as a script argument
if length(ARGS) == 1
    addprocs(parse(Int64, ARGS[1]))
elseif length(ARGS) > 1
    throw(ArgumentError("Only one command line argument (cores)."))
end

# Import packages needed for all workers
@everywhere  begin
    using Jevo
    using Distributions
    using DelimitedFiles
    using SharedArrays
    using LinearAlgebra
end

# Parameters
@everywhere begin
    gap = 10
    l_0 = 20
    fl = .7l_0
    f0 = 20l_0
    κ_arr = 0:2:20
    n = 4
    N = 1000
    steps = 10^8
    reps = 200
    F = Jevo.num_fermi(n, l_0, gap, f0/2N, fl/2N)
    emat = gap/l_0 * (ones(n, n) - Matrix{Float64}(I, n, n))
end

# Run one rep
@everywhere function run(N, l, emat, F, κ, l_0, gap, steps)
    pop = Jevo.mono_pop(N=N, l=l)
    Jevo.initiate!(pop, emat)
    for i in 1:steps
        Jevo.bp_substitution!(pop, emat, F)
        if rand() < κ/N
            Jevo.driver_mutation!(pop)
        end
        if (rand() < 1/10N) && (Jevo.get_energy(pop, emat)*l_0/length(pop.seqs)/gap < Jevo.γ_0(n))
            Jevo.l_substitution!(pop, emat, F)
	elseif (Jevo.get_energy(pop, emat)*l_0/length(pop.seqs)/gap > Jevo.γ_0(n))
            pop = Jevo.mono_pop(N=N, l=length(pop.seqs))
	    Jevo.initiate!(pop, emat)
	end
    end
    return Jevo.get_energy(pop, emat) * l_0/length(pop.seqs)/gap, length(pop.seqs)
end

# Store Metadata
open(date*"_METADATA.txt", "a") do io
    write(io, "gap=$gap\n")
    write(io, "l_0=$l_0\n")
    write(io, "f0=$f0\n")
    write(io, "fl=$fl\n")
    write(io, "kappa=$κ_arr\n")
    write(io, "n=$n\n")
    write(io, "N=$N\n")
    write(io, "steps=$steps\n")
    write(io, "reps=$reps")
end

E_results = SharedArray{Float64, 2}(length(κ_arr), reps)
l_list = SharedArray{Float64, 2}(length(κ_arr), reps)
kappa_list = SharedArray{Float64, 2}(length(κ_arr), reps)

# Run simulations and enjoy speed
@sync @distributed for j in 1:reps
    for (i1, κ) in enumerate(κ_arr)
        E, l= run(N, 150, emat, F, κ, l_0, gap, steps)
        E_results[i1, j] = E
        l_list[i1, j] = l
        kappa_list[i1, j] = κ
    end
    println("Run $j done.")
end
df = DataFrame(gamma=[(E_results...)...], l=[(l_list...)...], kappa=[(kappa_list...)...])
CSV.write(date * "_results.csv", df)


text:
module TestEnumerations

using  NFLTables
using  Test

@testset "enumeration tests" begin
    @test parse(SeasonPart, "PRE")  === PRE
    @test parse(SeasonPart, "REG")  === REG
    @test parse(SeasonPart, "POST") === POST

    @test length(instances(SeasonPart)) == 3
    @test length(instances(SuperBowl)) == 56

    @test Int(NFLTables.SB_LIV) == 2019
    @test SuperBowl(2019) === NFLTables.SB_LIV
end

end  # module


text:
using Pkg

Pkg.activate(@__DIR__)
Pkg.instantiate()

Pkg.status()

using PkgBenchmark

results = benchmarkpkg(
    dirname(@__DIR__),
    BenchmarkConfig(
        env = Dict(
            "JULIA_NUM_THREADS" => "1",
            "OMP_NUM_THREADS" => "1",
        ),
    )
)

dir_results = joinpath(@__DIR__, "results")
mkpath(dir_results)

writeresults(joinpath(dir_results, "$(results.date).json"), results)
export_markdown(joinpath(dir_results, "$(results.date).md"), results)


text:
# Analyse multisession solution progression

using Caesar
# IIF and KDE should not be necessary, but adding for late night run safety
using IncrementalInference
using KernelDensityEstimate

# get credentials
addrdict = getcredentials(nparticles=true, multisession=true)

# setup cloudgraphs interface
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)

# extras
multisessions = addrdict["multisession"]
N = parse(Int,addrdict["num particles"])

# RESET ENTIRE SYSTEM
removeMultisessions!(cloudGraph, session=multisessions[1])
removeMultisessions!(cloudGraph, session=multisessions[2])
removeMultisessions!(cloudGraph, session=multisessions[3])

addrdict["session"] = multisessions[1]
resetconvertdb(addrdict=addrdict)
setDBAllReady!(cloudGraph.neo4j.connection, addrdict["session"])
convertdb(addrdict=addrdict, N=N)
addrdict["session"] = multisessions[2]
resetconvertdb(addrdict=addrdict)
setDBAllReady!(cloudGraph.neo4j.connection, addrdict["session"])
convertdb(addrdict=addrdict, N=N)
addrdict["session"] = multisessions[3]
resetconvertdb(addrdict=addrdict)
setDBAllReady!(cloudGraph.neo4j.connection, addrdict["session"])
convertdb(addrdict=addrdict, N=N)

resultsfolder = "results_$(Dates.now(Dates.UTC))"
run(`mkdir -p $(resultsfolder)`)

for ANAITER in Int[1;2;3;4;5]
  for session in multisessions
    @show addrdict["session"] = session
    slamindb(addrdict=addrdict, iterations=1)
  end
  for session in multisessions
    multisessionsl = Vector{String}(setdiff(multisessions, [session]))
    # These are the base steps in rmInstMultisessionPriors!
    lm2others = getLandmOtherSessNeoIDs(cloudGraph,
                    session=session,multisessions=multisessionsl)
    # grab local subgraph using NeoIDs
    sfg, lms = getLocalSubGraphMultisession(cloudGraph, lm2others,
                    session=session, numneighbors=1)

    fg = initfg(sessionname=session, cloudgraph=cloudGraph)
    fullLocalGraphCopy!(fg)

    jldfile = string(resultsfolder,"/", session, "_$(ANAITER).jld")
    savejld(fg, file=jldfile)
    for sym in lms
      @show sym
      fv = getVert(fg, sym, nt=:fnc)
      ffv = getData(fv).fnc.usrfnc!
      # plotKDE(ffv.belief)

      p = getVertKDE(fg, sym)
      pl = plotKDE([ffv.belief; p], c=["red";"green"], levels=3, legend=["MultiSess"; "Current"], title=string(sym), fill=true)

      filename = string(resultsfolder,"/", session,"_",sym, "_$(ANAITER)")
      saveplot(pl, name=filename, frt=:pdf, nw=true)
      saveplot(pl, name=filename, frt=:png, nw=true)
    end
  end
end


text:
# provide a way to invoke any solver

"""
  This function provides a way to invoke any physics solver based on the
  specification of the physics in the input file.

  The physics module must have already been registered using [`register_physics`](@ref)

  This function is equivalent to [`solvePDE`](@ref), it is maintained here
  for legacy purposes only.

  **Inputs**

   * input_file: an AbstractString specifying the path to the input file

  **Outputs**

   * mesh: the AbstractMesh object used during the solve
   * sbp: the SBP operator used by the solver
   * eqn: the AbstractSolutionData object during the solve.  At exit,
          eqn.q_vec should have the final solution in it
   * opts: the options dictionary
"""
function run_solver(input_file::Union{AbstractString, Dict})

  return solvePDE(input_file)

end


text:
##
# These routines simplify the multiplication of two operators
# they are currently used for dekron of PDEs to make sure equivalent operators are
# recognized as the same
# I put this in Extras as it's currently a "hack"
#
# In the future, product rule/etc. could be implemented, and some sort of ordering
# to recognize M[f,C^(1)] Conversion(T,C^(1))==Conversion(T,C^(1))M[f,T]
##

simplifytimes(A::Operator,B::Operator) = isconstop(A)||isconstop(B)?A*B:[A;B]
simplifytimes(A::Operator,B::Vector) = [simplifytimes(A,B[1]);B[2:end]]
simplifytimes(A::Vector,B::Operator) = [A[1:end-1];simplifytimes(A[end],B)]
function simplifytimes(A::ConstantOperator,B::ConstantOperator)
    @assert A.space == B.space
    ConstantOperator(Number(A)*Number(B),A.space)
end
function simplifytimes(A::SpaceOperator,B::SpaceOperator)
   @assert domainspace(A)==rangespace(B)
    SpaceOperator(simplifytimes(A.op,B.op),domainspace(B),rangespace(A))
end

simplifytimes(A::Conversion,B::Conversion) = Conversion(domainspace(B),rangespace(A))

simplify(A)=A

function simplify(A::TimesOperator)
    ret=foldr(simplifytimes,A.ops)
    if isa(ret,Vector)
        @assert length(ret)>1
        TimesOperator(ret)
    else
        ret
    end
end


text:
easy_print(w::Ptr{Void}, y, x, str::ASCIIString) = TermWin.mvwprintw(w,y - 1,x - 1,"%s",str)

function print_cell(b::Board, w::Ptr{Void}, y, x)
    b.matrix_rep[b.active,y,x,1] == 1 && (easy_print(w, y+1, x+1, "@"); return)
    b.matrix_rep[b.active,y,x,2] == 1 && (easy_print(w, y+1, x+1, "+"); return)
    b.matrix_rep[b.active,y,x,3] == 1 && (easy_print(w, y+1, x+1, "*"); return)
end

function print_corners(b::Board,w::Ptr{Void})
    easy_print(w,1,1, "+")
    easy_print(w,1,b.width+2, "+")
    easy_print(w,b.height+2,1, "+")
    easy_print(w,b.height+2,b.width+2, "+")
end

function print_border(b::Board,w::Ptr{Void})
    for i=1:b.width
        easy_print(w,1,i+1,"-")
        easy_print(w,b.height+2,i+1,"-")
    end
    for i=1:b.height
        easy_print(w,i+1,1, "|")
        easy_print(w,i+1,b.width+2, "|")
    end
end

function print_frame(b::Board,w::Ptr{Void})
    print_corners(b,w)
    print_border(b,w)
end

function print_field(b::Board, w::Ptr{Void})
    for y = 1:b.height, x=1:b.width
        print_cell(b,w,y,x)
    end
end

max_stars_repo_name: RatanRSur/robots-from-scratch

max_issues_repo_name: RatanRSur/robots-from-scratch

max_forks_repo_name: RatanRSur/robots-from-scratch


text:
# audio.jl

@testset verbose= true "audio.jl" begin

    # Example audio
    audio_one_channel = SampleBuf(rand(1000), 10)
    audio_two_channel = SampleBuf(rand(1000,2), 10)
    audio_multi_channel = SampleBuf(rand(1000,2), 10)

    # Tests for mono
    @testset "mono" begin
        @test size(mono(audio_two_channel)) == (1000,)
        @test size(mono(audio_multi_channel)) == (1000,)

    end

    # Tests for resample method
    @testset "resample.jl" begin
        @test size(resample(audio_one_channel, 5))[1] == 500
        @test size(resample(audio_two_channel, 5))[1] == 500
        @test size(resample(audio_multi_channel, 5))[1] == 500
    end

    # Tests for duration method
    @testset "duration" begin
        @test duration(audio_one_channel) == 100
        @test duration(audio_two_channel) == 100
        @test duration(audio_multi_channel) == 100
    end
    audio_one_channel = SampleBuf(rand(10000), 1000)
    audio_two_channel = SampleBuf(rand(10000,2), 1000)
    audio_multi_channel = SampleBuf(rand(10000,2), 1000)
    # Tests for pitchshift method
    @testset "pitchshift" begin
        @test isapprox(duration(pitchshift(audio_one_channel)), duration(audio_one_channel), atol=0.2)
        @test isapprox(duration(pitchshift(audio_two_channel)), duration(audio_two_channel), atol=0.2)
        @test isapprox(duration(pitchshift(audio_multi_channel)), duration(audio_multi_channel), atol=0.2)
    end

    audio_one_channel = SampleBuf(rand(10000), 1000)
    audio_two_channel = SampleBuf(rand(10000,2), 1000)
    audio_multi_channel = SampleBuf(rand(10000,2), 1000)
    # Tests for speedup method
    @testset "speedup" begin

        @test isapprox(duration(speedup(audio_one_channel, 2)), duration(audio_one_channel) / 2.0 , atol=0.4)
        @test isapprox(duration(speedup(audio_two_channel, 2)), duration(audio_two_channel) / 2.0 , atol=0.4)
        @test isapprox(duration(speedup(audio_multi_channel, 2)), duration(audio_multi_channel) / 2.0 , atol=0.4)
    end

    audio_one_channel = SampleBuf(rand(10000), 1000)
    audio_two_channel = SampleBuf(rand(10000,2), 1000)
    audio_multi_channel = SampleBuf(rand(10000,2), 1000)
    # Tests for slowdown method
    @testset "slowdown" begin

        @test isapprox(duration(slowdown(audio_one_channel, 2)), duration(audio_one_channel) * 2.0, atol=1.2)
        @test isapprox(duration(slowdown(audio_two_channel, 2)), duration(audio_two_channel) * 2.0, atol=1.2)
        @test isapprox(duration(slowdown(audio_multi_channel, 2)), duration(audio_multi_channel) * 2.0, atol=1.2)
   end

    # Tests for zero crossing rate method
    @testset "zero_crossing_rate" begin
        audio_one_channel = SampleBuf(rand(100000), 10)

        expected = length(zero_crossing_rate(audio_one_channel, 2048, 512))
        res = 192
        @test expected == res
    end

end


text:
using FiniteDifferences
using JET
using Lux
using Random
using Zygote

function Base.isapprox(nt1::NamedTuple{fields}, nt2::NamedTuple{fields};
                       kwargs...) where {fields}
    checkapprox(xy) = isapprox(xy[1], xy[2]; kwargs...)
    checkapprox(t::Tuple{Nothing, Nothing}) = true
    return all(checkapprox, zip(values(nt1), values(nt2)))
end

function Base.isapprox(t1::NTuple{N, T}, t2::NTuple{N, T}; kwargs...) where {N, T}
    checkapprox(xy) = isapprox(xy[1], xy[2]; kwargs...)
    checkapprox(t::Tuple{Nothing, Nothing}) = true
    return all(checkapprox, zip(t1, t2))
end

# Test the gradients generated using AD against the gradients generated using Finite Differences
function test_gradient_correctness_fdm(f::Function, args...; kwargs...)
    gs_ad = Zygote.gradient(f, args...)
    gs_fdm = FiniteDifferences.grad(FiniteDifferences.central_fdm(5, 1), f, args...)
    for (g_ad, g_fdm) in zip(gs_ad, gs_fdm)
        @test isapprox(g_ad, g_fdm; kwargs...)
    end
end

# Some Helper Functions
function run_fwd_and_bwd(model, input, ps, st)
    y, pb = Zygote.pullback(p -> model(input, p, st)[1], ps)
    gs = pb(ones(eltype(y), size(y)))
    # if we make it to here with no error, success!
    return true
end

function run_model(m::Lux.AbstractExplicitLayer, x, mode=:test)
    ps, st = Lux.setup(Random.default_rng(), m)
    if mode == :test
        st = Lux.testmode(st)
    end
    return Lux.apply(m, x, ps, st)[1]
end

# JET Tests
function run_JET_tests(f, args...; call_broken=false, opt_broken=false, kwargs...)
    @static if VERSION >= v"1.7"
        test_call(f, typeof.(args); broken=call_broken)
        test_opt(f, typeof.(args); broken=opt_broken, target_modules=(Lux,))
    end
end


text:
# ===========================================================================
# Helicopter
#
# system type: linear continuous system
# state dimension: 28
#
# This is an 8-dimensional helicopter model with a 20-dimensional controller.
#
# [1] S. Skogestad and I. Postlethwaite: Multivariable Feedback Control:
# Analysis and Design. John Wiley & Sons, 2005.
# ===========================================================================
using MathematicalSystems, SparseArrays

function helicopter_A()
    A = [[0, 0, 0, 0.998574, 0.0533843, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, -0.00318222, 0.0595247, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, -11.5705, -2.54464, -0.0636026, 0.106781, -0.0949187, 0.00710757, -5.20818, -23.1259, -4.24693, -0.710744, -1.4867, -0.0445482, 0.0343664, 0.0229902, 0.580153, 1.04755, 11.265, 2.62228, 0.382949, 0.453076, 0.0226519, 0.0125336, 0.108399, -0.160139, -6.11536, -0.0389724],
    [0, 0, 0.439357, -1.99818, 0, 0.0166519, 0.018462, -0.00118747, -10.7788, 0.343053, -0.0349493, -3.46674, 0.0134159, -0.00650175, -0.00724803, -0.00483591, -0.152495, 2.03915, -0.169614, -0.279125, -0.00367822, -0.00167971, 0.050812, 0.0432735, -0.0623053, 0.979163, 0.0404014, 0.010697],
    [0, 0, -2.0409, -0.458999, -0.735028, 0.0192558, -0.00459562, 0.00212036, -1.06433, -3.65916, -0.493547, -0.0595611, -1.95351, -0.0109391, -0.000271004, 0.0505639, 1.58834, 0.219924, 1.82192, 2.69319, 0.0686196, 0.08631, 0.00241543, 0.000727041, 0.640964, -0.0230167, -1.06985, -0.596675],
    [-32.1036, 0, -0.503355, 2.29786, 0, -0.0212158, -0.0211679, 0.0158116, 12.4592, -0.353637, 0.0641365, 4.03533, -0.272153, 0.00099959, 0.00252969, 0.111633, 2.73642, -2.35454, 0.175217, 0.863082, 0.00440128, 0.00437341, -0.059795, -0.0510055, 0.609532, -1.11578, -0.0263051, -0.0337389],
    [0.102161, 32.0578, -2.34722, -0.503612, 0.834948, 0.0212266, -0.0378797, 0.000354004, -0.76148, -5.12933, -1.13138, -0.148788, 1.38074, -0.012931, 0.00764572, 0.0737967, 1.36175, 0.150701, 2.45244, -1.44884, 0.0765163, 0.0871226, 0.00462319, 0.00263523, -0.0794019, -0.0310234, -1.22553, 0.657926],
    [-1.91097, 1.71383, -0.00400543, -0.0574112, 0, 0.0139896, -0.000906753, -0.290514, -2.06064, -0.774846, -0.471752, -1.21389, 5.03052, 0.126408, 0.113189, -2.07842, -50.1852, 0.340666, 0.375863, -10.6412, -0.00363415, -0.0479628, 0.0305097, 0.0275842, -10.5424, -0.126841, -0.391839, 0.420788],
    [0.126296, -0.00289825, -0.319316, 0.785202, 0.00177237, 5.84372e-06, 5.23381e-06, -9.78995e-05, -0.0726115, 0.00166629, 0.195701, 0.517339, 0.0523653, -3.35973e-06, -3.00908e-06, 5.62854e-05, -0, -0, -0, -0, -0, 0, -0, 0, 0, 0, 0, 0],
    [-0.0181141, 0.077615, 0.710897, 2.45428, -0.012793, 4.05106e-05, 3.62825e-05, -0.000678672, 0.0104143, -0.0446232, 0.564308, -1.50732, 0.0668797, -2.32908e-05, -2.08599e-05, 0.000390189, 0, 0, -0, -0, -0, 0, -0, 0, 0, 0, 0, 0],
    [-0.0199573, 0.00710897, 122.639, 1.7917, 0.13833, 7.26169e-07, 6.50379e-07, -1.21655e-05, -8.48115, -37.7139, -93.6582, -4.80197, -2.56739, 0.0341383, -0.0388801, 0.0446032, 0.946017, 1.70817, 18.3691, 4.27597, 0.62445, 0.738801, 0.0369369, 0.0204377, 0.17676, -0.261129, -9.9719, -0.0635496],
    [0.00785296, 0.00392575, 0.287856, 58.0535, 0.0306981, -6.8376e-06, -6.12396e-06, 0.00011455, -17.5807, 0.557136, 0.205947, -43.2308, 0.00422708, 0.00605385, 0.00664669, -0.00913893, -0.248663, 3.32511, -0.276579, -0.455151, -0.00599782, -0.00273899, 0.0828557, 0.0705632, -0.101597, 1.59665, 0.0658798, 0.0174429],
    [0.0114973, -0.012583, 13.8484, 22.2888, 0.0422873, 0.000197559, 0.00017694, -0.00330969, -1.74214, -5.95951, -11.333, -14.2165, -3.9448, 0.00130458, -0.00513926, 0.0864743, 2.59, 0.358615, 2.97089, 4.3916, 0.111893, 0.14074, 0.00393867, 0.00118554, 1.04518, -0.0375318, -1.74453, -0.972958],
    [-16.9391, 0.618054, 107.921, 204.524, 0.204408, 0.0047421, 0.00424717, -0.0794442, -2.04846, -0.93199, -66.5409, -116.47, -0.561301, -0.0223122, -0.0194847, 0.243519, 4.4621, -3.83939, 0.285714, 1.40737, 0.00717686, 0.00713142, -0.0975037, -0.0831712, 0.993922, -1.81943, -0.0428939, -0.0550157],
    [-0.54281, 7.08182, -135.012, 460.929, 0.0364986, 0.00693724, 0.0062132, -0.116219, -0.827455, 19.6222, 78.5537, -283.239, 3.06544, -0.00384762, -0.0289845, 0.187507, 2.22051, 0.245738, 3.99902, -2.36252, 0.12477, 0.142065, 0.00753873, 0.0042971, -0.129475, -0.0505877, -1.99839, 1.07284],
    [-1.28646, 0.142279, -1.26575, 65.7431, -1.3207, -0.061856, -0.0554001, 1.03627, -4.53149, 0.368539, 0.0024871, -42.3265, 8.96223, 0.255677, 0.215513, -4.27544, -81.8337, 0.5555, 0.612895, -17.3518, -0.00592597, -0.0782097, 0.0497501, 0.0449796, -17.1907, -0.206831, -0.638946, 0.686151],
    [0, 0, 0, 0, 0, -0.00970226, -0.00868964, 0.162541, 0, 0, 0, 0, 0, 0, 0, 0, -0.012, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [-8.15316, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.005, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, -3.26127, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.005, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0.174412, -3.26127, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.01, 0, 0, 0, 0, 0, 0, 0, 0],
    [-0, 0, -3.26127, -0, -0, -0, -0, 0, 0, -0, -0, 0, 0, 0, 0, -0, 0, 0, 0, -0, -8.5, -18, 0, 0, 0, 0, 0, 0],
    [0, 0, -0, -0, 0, -0, -0, 0, -0, 0, 0, 0, 0, 0, 0, -0, -0, -0, -0, -0, 1, 0, 0, 0, 0, 0, 0, 0],
    [-0, 0, 0, -8.15316, -0, -0, -0, 0, 0, 0, -0, 0, 0, 0, 0, -0, -0, 0, -0, 0, -0, 0, -8.5, -18, 0, 0, 0, 0],
    [0, 0, -0, -0, 0, -0, -0, 0, -0, 0, 0, 0, 0, 0, 0, -0, 0, -0, -0, -0, -0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0.699961, 0.262038, 0.15959, 0.411552, -1.70162, -0.0427567, -0.0382852, 0.703046, 16.9757, -0.115788, -0.127109, 3.59954, 0.00122974, 0.0162237, -0.010334, -0.00934236, -6.43393, 0.0426396, 0.132541, -0.142338],
    [0, 0, 0, 0, 0, 0, 0, 0, -37.0015, 0.783589, -0.183855, -11.8696, -0.106084, -0.0263066, -0.0271181, 0.036745, 0.764602, 7.00237, -0.390319, -0.642631, -0.00570167, 0.00352225, 0.173868, 0.147911, 0.0560927, -6.64198, 0.0396022, 0.0261817],
    [0, 0, 0, 0, 0, 0, 0, 0, 1.9914, 13.76, 2.53042, 0.0825288, 0.728265, 0.0239028, -0.0228966, 0.0153276, 0.370477, -0.412566, -6.70095, -1.32704, -0.227019, -0.267482, -0.00865099, -0.00339436, 0.0987926, 0.197714, -6.3694, -0.0119768],
    [0, 0, 0, 0, 0, 0, 0, 0, 1.96586, -3.74313, -1.96265, 0.0189294, 11.145, -0.036002, -0.00122215, 0.602488, 11.6398, -0.407673, 1.50774, -12.8, 0.0053931, -0.0142088, -0.000915158, -0.000640326, -0.0365353, 0.012459, -0.083125, -5.47283]]

    return sparse(hcat(A...)')
end

function helicopter_model()
    # system matrix
    A = helicopter_A()

    # linear continuous system
    S = @system(x' = Ax)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
import Random: AbstractRNG

# No info
"""
    Flat <: ContinuousUnivariateDistribution

A distribution with support and density of one
everywhere.
"""
struct Flat <: ContinuousUnivariateDistribution end

Distributions.rand(rng::AbstractRNG, d::Flat) = rand(rng)
Distributions.logpdf(d::Flat, x::Real) = zero(x)
Distributions.minimum(d::Flat) = -Inf
Distributions.maximum(d::Flat) = +Inf

# For vec support
Distributions.logpdf(d::Flat, x::AbstractVector{<:Real}) = zero(x)

"""
    FlatPos(l::Real)

A distribution with a lower bound of `l` and a density
of one at every `x` above `l`.
"""
struct FlatPos{T<:Real} <: ContinuousUnivariateDistribution
    l::T
end

Distributions.rand(rng::AbstractRNG, d::FlatPos) = rand(rng) + d.l
Distributions.logpdf(d::FlatPos, x::Real) = x <= d.l ? -Inf : zero(x)
Distributions.minimum(d::FlatPos) = d.l
Distributions.maximum(d::FlatPos) = Inf

# For vec support
function Distributions.logpdf(d::FlatPos, x::AbstractVector{<:Real})
    return any(x .<= d.l) ? -Inf : zero(x)
end

"""
    BinomialLogit(n<:Real, I<:Integer)

A univariate binomial logit distribution.
"""
struct BinomialLogit{T<:Real, I<:Integer} <: DiscreteUnivariateDistribution
    n::I
    logitp::T
end

"""
    BinomialLogit(n<:Real, I<:Integer)

A multivariate binomial logit distribution.
"""
struct VecBinomialLogit{T<:Real, I<:Integer} <: DiscreteUnivariateDistribution
    n::Vector{I}
    logitp::Vector{T}
end

function logpdf_binomial_logit(n, logitp, k)
    logcomb = -StatsFuns.log1p(n) - SpecialFunctions.lbeta(n - k + 1, k + 1)
    return logcomb + k * logitp - n * StatsFuns.log1pexp(logitp)
end

function Distributions.logpdf(d::BinomialLogit{<:Real}, k::Int)
    return logpdf_binomial_logit(d.n, d.logitp, k)
end

function Distributions.logpdf(d::VecBinomialLogit{<:Real}, ks::Vector{<:Integer})
    return sum(logpdf_binomial_logit.(d.n, d.logitp, ks))
end

"""
    OrderedLogistic(η::Any, cutpoints<:AbstractVector)

An ordered logistic distribution.
"""
struct OrderedLogistic{T1, T2<:AbstractVector} <: DiscreteUnivariateDistribution
   η::T1
   cutpoints::T2

   function OrderedLogistic(η, cutpoints)
        if !issorted(cutpoints)
            error("cutpoints are not sorted")
        end
        return new{typeof(η), typeof(cutpoints)}(η, cutpoints)
   end

end

function Distributions.logpdf(d::OrderedLogistic, k::Int)
    K = length(d.cutpoints)+1
    c =  d.cutpoints

    if k==1
        logp= -softplus(-(c[k]-d.η))  #logp= log(logistic(c[k]-d.η))
    elseif k<K
        logp= log(logistic(c[k]-d.η) - logistic(c[k-1]-d.η))
    else
        logp= -softplus(c[k-1]-d.η)  #logp= log(1-logistic(c[k-1]-d.η))
    end

    return logp
end

function Distributions.rand(rng::AbstractRNG, d::OrderedLogistic)
    cutpoints = d.cutpoints
    η = d.η

    K = length(cutpoints)+1
    c = vcat(-Inf, cutpoints, Inf)

    ps = [logistic(η - i[1]) - logistic(η - i[2]) for i in zip(c[1:(end-1)],c[2:end])]

    k = rand(rng, Categorical(ps))

    if all(ps.>0)
        return(k)
    else
        return(-Inf)
    end
end


text:
function demcz_sample(logobj, Zmat, opts::DEMCopt; prevrun=nothing)
    return demcz_sample(logobj, Zmat, opts.N, opts.K, opts.Ngeneration, opts.Nblocks, opts.blockindex, opts.eps_scale, opts.γ; prevrun=prevrun, verbose = opts.verbose, print_step=opts.print_step, autostop=opts.autostop, autostop_Rhat=opts.autostop_Rhat, autostop_every = opts.autostop_every)
end
"""
demcz_sample(logobj, Zmat, N=4, K=10, Ngeneration=5000, Nblocks=1, blockindex=[1:size(Zmat,2)], eps_scale=1e-4*ones(size(Zmat,2)), γ=2.38; prevrun=nothing, verbose = true, print_step=100)

Serial run of DEMC chain. At each generation all chains are updated one after the other.
"""
function demcz_sample(logobj, Zmat, N=4, K=10, Ngeneration=5000, Nblocks=1, blockindex=[1:size(Zmat,2)], eps_scale=1e-4*ones(size(Zmat,2)), γ=2.38; prevrun=nothing, verbose = true, print_step=100, autostop=:no, autostop_Rhat=1.01, autostop_every = 1000)
    nrowZ, d = size(Zmat)
    Zmat = vcat(Zmat, zeros(Int(ceil(N*Ngeneration/K)), d))
    M = ones(Int64,1)*nrowZ
    if prevrun == nothing
        # start chain at last N elements of Z
        X = Zmat[end-N+1:end, :]
        # initial obj values
        log_objcurrent = map(logobj, [X[i,:] for i = 1:N])
    else
        # start chain at last generation of prevrun
        X = prevrun.chain[:,:,end]
        log_objcurrent = prevrun.log_objcurrent[:,end]
    end
    # preallocate the chain
    mc = MC(zeros(N,  d, Ngeneration),zeros(N, Ngeneration), X, log_objcurrent)
    # print inital values
    if verbose
        print_status(mc, 0)
    end
    # run through generations
    for ig = 1:Ngeneration
        for ic = 1:N
            runchain!(ic, ig, ig, mc, Zmat, K, M, logobj, blockindex, eps_scale, γ, Nblocks)
        end
        if verbose
            if mod(ig, print_step) == 0.
                print_status(mc, ig;printlast=autostop_every)
            end
        end
        if autostop == :Rhat
            if mod(ig, autostop_every) == 0
                Rhat = Rhat_gelman(mc.chain[:,:, ig-autostop_every+1:ig], N, autostop_every, d)
                accept_ratio = sum(diff(mc.log_obj[ig-autostop_every+1:ig], dims = 2).!=0., dims = 2)./ (autostop_every-1)
                if (maximum(Rhat) < autostop_Rhat)
                    if mean(accept_ratio) < 0.1
                        println("Warning: accept ratio below 10% on average")
                    end
                    mc = MC(mc.chain[:,:,1:ig],  mc.log_obj[:, 1:ig], mc.Xcurrent, mc.log_objcurrent)
                    if prevrun != nothing
                        return MC(cat(prevrun.chain, mc.chain, dims=3),cat(prevrun.log_obj, mc.log_obj, dims=2), mc.Xcurrent, mc.log_objcurrent), Zmat
                    else
                        return mc, Zmat[1:M[1],:]
                    end
                end
            end
        end
    end

    if prevrun != nothing
        return MC(cat(prevrun.chain, mc.chain, dims=3),cat(prevrun.log_obj, mc.log_obj, dims=2), mc.Xcurrent, mc.log_objcurrent), Zmat
    else
        return mc, Zmat
    end
end

function print_status(mc, ig;printlast=500)
    #bestval = maximum(mc.log_objcurrent)
    #bestpar = mc.Xcurrent[findfirst(bestval.==mc.log_objcurrent), :]
    avgpar = mean(mc.chain[:,:,max(1,ig-printlast):max(1,ig)], dims = [1,3])
    avgval = mean(mc.log_obj[:,max(1,ig-printlast):max(1,ig)], dims = [1,2])

    println("-----------------------")
    println("iteration $ig")
    println("average par = $avgpar")
    println("average val = $avgval")
    #println("bestval = $bestval")
    #println("bestpar = $bestpar")
    println("-----------------------")
end

function runchain!(ic, from, to, mc, Zmat, K, M, logobj, blockindex, eps_scale, γ, Nblocks)
    for ig = from:to
        Xcurrent, current_logobj = update_blocks(mc.Xcurrent[ic, :], mc.log_objcurrent[ic], Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
        # update in chain
        mc.chain[ic, :, ig] .= Xcurrent
        mc.log_obj[ic, ig] = current_logobj
        mc.Xcurrent[ic, :] .= Xcurrent
        mc.log_objcurrent[ic] = current_logobj
        if mod(ig, K) == 0.
            Zmat[M[1]+1,:] .= mc.Xcurrent[ic,:]
            M .= M .+ 1
        end
    end
end

"""
demcz_sample_par(logobj, Zmat, opts::DEMCopt; sync_every = 1000, prevrun=nothing)

Runs each chain on a separate process - Z is updated simultaenously among all chains running in parallel.
"""
function demcz_sample_par(logobj, Zmat, opts::DEMCopt; sync_every = 1000, prevrun=nothing, workerpool = CachingPool(workers()))
    return demcz_sample_par(logobj, Zmat, opts.N, opts.K, opts.Ngeneration, opts.Nblocks, opts.blockindex, opts.eps_scale, opts.γ;sync_every = sync_every, prevrun=prevrun,  autostop=opts.autostop, autostop_Rhat=opts.autostop_Rhat, verbose=opts.verbose, workerpool=workerpool)
end

function demcz_sample_par(logobj, Zmat, N=4, K=10, Ngeneration=5000, Nblocks=1, blockindex=[1:size(Zmat,2)], eps_scale=1e-4*ones(size(Zmat,2)), γ=2.38; sync_every = 1000, prevrun=nothing, autostop=:Rhat, autostop_Rhat=1.1, verbose=true, workerpool = CachingPool(workers()))
    # prep storage etc
    nrowZ, d = size(Zmat)
    global Zshared = SharedArray(vcat(Zmat, zeros(Int(ceil(N*Ngeneration/K)), d)))
    global M = SharedArray(ones(Int64,1)*nrowZ)
    if prevrun == nothing
        # start chain at last N elements of Z
        X = Zshared[M[1]-N+1:M[1], :]
        # initial obj values
        log_objcurrent = pmap(logobj, workerpool, [X[i,:] for i = 1:N])
    else
        # start chain at last generation of prevrun
        X = prevrun.chain[:,:,end]
        log_objcurrent = prevrun.log_objcurrent[:,end]
    end

    # preallocate the chain
    global mc = MCShared(SharedArray(zeros(N,  d, Ngeneration)),SharedArray(zeros(N, Ngeneration)), SharedArray(X), SharedArray(log_objcurrent))
    global Ngenerationg = Ngeneration
    # define on all cores
    passobj(myid(), workers(), [:Zshared,:M, :Ngenerationg], from_mod=DEMC, to_mod=DEMC)
    passobj(myid(), workers(), [:mc], from_mod=DEMC, to_mod=DEMC)
    # split up generations
    Nsets = ceil(Ngeneration/sync_every)
    splitgens = [(Int((i-1)*sync_every+1),Int(min(i*sync_every, Ngeneration)))  for i = 1:Nsets]
    global from = 0
    global to = 0
    for set in splitgens
        global from = set[1]
        global to = set[2]
        passobj(myid(), workers(), [:from, :to], from_mod=DEMC, to_mod=DEMC)
        pmap(ic -> runchain!(ic, from, to, mc, Zshared, K, M, logobj, blockindex, eps_scale, γ, Nblocks), workerpool, 1:N)
        if verbose
            print_status(mc, to; printlast=sync_every)
        end
        if autostop == :Rhat
            Rhat = Rhat_gelman(mc.chain[:,:, from:to], N, to-from, d)
            accept_ratio = sum(diff(mc.log_obj, dims = 2).!=0., dims = 2)./ (sync_every-1)
            if (maximum(Rhat) < autostop_Rhat)
                if mean(accept_ratio) < 0.1
                    println("Warning: accept ratio below 10% on average")
                end
                mc = MC(mc.chain[:,:,1:to],  mc.log_obj[:, 1:to], mc.Xcurrent, mc.log_objcurrent)
                if prevrun != nothing
                    mc = MCShared(cat(prevrun.chain, mc.chain, dims=3),cat(prevrun.log_obj, mc.log_obj, dims=2), mc.Xcurrent, mc.log_objcurrent)
                    return mc, Zshared[1:M[1],:]
                else
                    return mc, Zshared[1:M[1],:]
                end
            end
        end
    end

    if prevrun != nothing
        mc = MCShared(cat(prevrun.chain, mc.chain, dims=3),cat(prevrun.log_obj, mc.log_obj, dims=2), mc.Xcurrent, mc.log_objcurrent)
        return mc, Zmat
    else
        return mc, Zmat
    end
end

function update_blocks(Xcurrent, current_logobj, Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
    for ib in 1:Nblocks
        Xcurrent, current_logobj = update_demcz_chain_block(Xcurrent, current_logobj, ib, Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
    end
    return Xcurrent, current_logobj
end

function update_demcz_chain_block(Xcurrent, current_logobj, ib, Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
    # generate proposal
    set = collect(1:M[1])
    i1 = rand(set)
    deleteat!(set, i1)
    i2 = rand(set)
    de_diffvec = zeros(eltype(Xcurrent), length(Xcurrent))
    block = blockindex[ib]
    blocklen = length(block)
    if blocklen == 1
        de_diffvec[block] = γ*(Zmat[i1,block]- Zmat[i2,block]) .+ eps_scale[block] .* randn()
    else
        de_diffvec[block] = γ/sqrt(2*blocklen)*(Zmat[i1,block]-Zmat[i2,block]) .+ eps_scale[block] .* randn(blocklen)
    end
    Xproposal = Xcurrent .+ de_diffvec
    log_objXprop = logobj(Xproposal)
    if accept(log_objXprop, current_logobj)
        return Xproposal, log_objXprop
    else
        return Xcurrent, current_logobj
    end
end

function accept(objprop, objprev)
    if log(rand()) < objprop - objprev
        return true
    else
        return false
    end
end


text:
module EffectsArmswingGaitAsymmetryStability

using Biomechanics, MAT, ChaosTools, Interpolations, UnsafeArrays, DSP

# Standard library
using Statistics, DelimitedFiles, LinearAlgebra

export readsegment,
       readtrials,
       analyzetrial

export DSSteadyState,
       SteadyStateSeg

abstract type DSSteadyState <: AbstractDataSource end

struct SteadyStateSeg <: DSSteadyState
    events::Dict{Symbol,Vector{Float64}}
    data::Dict{Symbol,Matrix{Float64}}
end

SteadyStateSeg() = SteadyStateSeg(Dict{Symbol,Matrix}(), Dict{Symbol,AbstractVector}())

function readsegment(DSData::Type{<:AbstractDataSource},
                     trial::Trial{<:AbstractDataSource},
                     st::Float64 = 0.0,
                     en::Float64 = Inf;
                     events::Vector = [],
                     ts::Vector = [],
                     fs = 100)
    isfile(trial.path) || throw(ArgumentError("trial $(trial.path) exist"))
    st >= 0.0 || throw(ArgumentError("start time must be positive"))
    st <= en || throw(ArgumentError("end time must be greater than start time"))

    evnames = Symbol.(events)
    file = matopen(trial.path)

    revents = Dict{Symbol,Vector{Float64}}()
    for e in events
        if exists(file, string(e))
            syme = Symbol(e)
            tmp = read(file, string(e))[1]
            if tmp isa AbstractArray
                revents[syme] = vec(tmp)
            else
                revents[syme] = [tmp]
            end
            strt = findfirst(x -> x >= st, revents[syme])
            if strt == 0
                # @warn "no $e events during given start and end times"
                delete!(revents, syme)
                break
            end
            endi = findlast(x -> x <= en, revents[syme])
            revents[syme] = revents[syme][strt:endi] .- st .+ (1/fs) # Shift events to be index accurate for data subsection
        else
            # @warn "Requested event $e does not exist in source data"
        end
    end

    data = Dict{Symbol,Matrix{Float64}}()
    for t in ts
        if exists(file, string(t))
            symt = Symbol(t)
            data[symt] = read(file, string(t))[1]
            len = size(data[symt], 1)
            strti = round(Int, st*fs)
            endi = en == Inf ? len : min(len, round(Int, en*fs))
            data[symt] = data[symt][strti:endi, :]
        else
            @warn "Requested time series $t does not exist in source data"
        end
    end

    return Segment(trial, Dict{Symbol,Any}(), DSData(revents, data))
end

const fs = 100

function readtrials(datadir::String)
    sstrials = Vector{Trial}()

    fp = abspath(datadir)

    trials = readdir(datadir)
    filter!(trial -> endswith(trial, "mat"), trials)

    for trial in trials
        conds = Dict{Symbol, Symbol}()
        m = match(r"S(?<sub>\d{1,2})_(?<arms>(held|norm|active))_(?<sym>(sym|asym))", trial)
        sub = tryparse(Int, m[:sub])
        conds[:sym] = Symbol(m[:sym])
        conds[:arms] = Symbol(m[:arms])

        push!(sstrials, Trial{DSSteadyState}(sub, splitext(basename(trial))[1], joinpath(fp, trial), conds))
    end

    return sstrials
end

"""
    continuousphase(θ, events) -> Vector

Calculate the continuous phase of `θ` by finding the angle of the Hilbert transform of the
amplitude centered signal, `θ`. The amplitude centering uses the average extrema between
pairs of `events`.

This method of continuous phase calculation is as recommended by Lamb and Stöckl (2014) [1]

[1] P. F. Lamb and M. Stöckl, “On the use of continuous relative phase: Review of current
    approaches and outline for a new standard,” Clinical Biomechanics, vol. 29, no. 5,
    pp. 484–493, May 2014, [doi](https://doi.org/10.1016%2Fj.clinbiomech.2014.03.008).
"""
function continuousphase(θ, events::AbstractVector{<:Integer})
    mi, ma = avgextrema(θ, events)
    θcent = θ .- Ref(mi) .- Ref((ma - mi)/2)
    Hθ = hilbert(θcent)

    return angle.(Hθ)
end

"""
    meanrange(x, events) -> Tuple(avg, std)

Find the average range and range variability of `x` for the set of all intervals given by
`events`.
"""
function meanrange(x::AbstractVector{T}, events::AbstractVector{<:Integer}) where {T}
    l = length(events) - 1
    E = Array{T}(undef, l, 2)
    _rangedextrema!(E, x, events, l)

    roms = diff(E; dims=2)

    avg = mean(roms)
    sd = std(roms)

    (avg, sd)
end

"""
    avgextrema(x, events) -> Tuple(min, max)

Find the average minima and maxima for the set of all intervals given by `events`.
"""
function avgextrema(x::AbstractVector{T}, events::AbstractVector{<:Integer}) where {T}
    l = length(events) - 1
    E = Array{T}(undef, l, 2)
    _rangedextrema!(E, x, events, l)

    avgex = mean(E; dims=1)

    (avgex[1], avgex[2])
end

function _rangedextrema!(E, x, events, l)
    @inbounds for i in 1:l
        GC.@preserve E x events i begin
            E[i,1], E[i,2] = extrema(uview(x, events[i]:(events[i+1]-1)))
        end
    end
    return E
end

# Asymmetry function (see Plotnik et al. 2005 doi:10.1002/ana.20452)
"""
    asymmetry(l,r)

Calculate the asymmetry between `l` and `r` using the equation used by Plotnik et al. (2005)
[1].

[1] M. Plotnik, N. Giladi, Y. Balash, C. Peretz, and J. M. Hausdorff, “Is freezing of gait
    in Parkinson’s disease related to asymmetric motor function?,” Ann Neurol., vol. 57,
    no. 5, pp. 656–663, May 2005 [doi](https://doi.org/10.1002%2Fana.20452).
"""
asymmetry(l,r) = abs(log(l/r))*100

"""
    circmean(x)

Compute the circular mean of `x` in radians.

[1] N. I. Fisher, Statistical Analysis of Circular Data. Cambridge University Press, 1993.
"""
function circmean(x)
    s = mean(sin, x)
    c = mean(cos, x)

    atan(s, c)
end

"""
    circstd(x)

Compute the circular standard deviation of `x` in radians.

[1] N. I. Fisher, Statistical Analysis of Circular Data. Cambridge University Press, 1993.
"""
function circstd(x)
    s = mean(sin, x)
    c = mean(cos, x)

    sqrt(-2*log(hypot(c, s)))
end

function analyzetrial(trial::Trial, numstrides::Int)
    cols = [ "TrunkLinVel", "TrunkAngVel", "RFootPos", "LFootPos", "LHip", "RHip", "LShoulder", "RShoulder" ]
    seg = readsegment(SteadyStateSeg, trial, 25.0; events=[:RFC, :LFC], ts=cols)
    insuffstrerr = "insufficient number of strides in trial $(trial.name), subject $(trial.subject)"
    length(seg.data.events[:RFC]) < numstrides+1 && throw(ArgumentError(insuffstrerr))

    results = Dict{Symbol,Any}()

    ########################################
    ## Interleaving steps
    ########################################

    # Gather rfc and lfc info into a NamedTuple, with named fields `time`, `pos`, and `leg`, (`:RFC` or `:LFC`)
    rfc = [ (time=seg.data.events[:RFC][i],
             pos=SVector{2}(seg.data.data[:RFootPos][
                            round(Int, seg.data.events[:RFC][i]*fs)
                            ,1:2]),
             leg=:RFC) for i in 1:(numstrides+1) ]
    if seg.data.events[:LFC][1] < seg.data.events[:RFC][1]
        lfc = [ (time=seg.data.events[:LFC][i],
                 pos=SVector{2}(seg.data.data[:LFootPos][
                                round(Int, seg.data.events[:LFC][i]*fs)
                                ,1:2]),
                 leg=:LFC) for i in 2:(numstrides+2) ]
    else
        lfc = [ (time=seg.data.events[:LFC][i],
                 pos=SVector{2}(seg.data.data[:LFootPos][
                                round(Int, seg.data.events[:LFC][i]*fs)
                                ,1:2]),
                 leg=:LFC) for i in 1:(numstrides+1) ]
    end

    # Sort all steps and check that they are alternating--no double steps--and that they are
    # in the expected locations:
    # Odd elements need to be RFC for later assumptions made on whether the odd/even elements
    # of the diff are right or left steps
    steps = sort([rfc; lfc], by=(x -> x.time))

    while any(x -> x.leg === :LFC, steps[1:2:end])
        badstep = min(something(findfirst(x -> x.leg === :LFC, steps[1:2:end]), typemax(Int)÷2)*2-1,
                      something(findfirst(x -> x.leg === :RFC, steps[2:2:end]), typemax(Int)÷2)*2)
        steps = deleteat!(steps, badstep)
    end
    all(x -> x.leg === :RFC, steps[1:2:end]) || throw(DomainError("Odd elements aren't all RFC's"))
    all(x -> x.leg === :LFC, steps[2:2:end]) || throw(DomainError("Even elements aren't all LFC's"))

    ########################################
    ## Gait variability
    ########################################

    # After `diff`ing, odd elements will be left steps, even elements will be right steps
    steptimes = diff(getindex.(steps, :time))
    if trial.conds[:sym] == :asym
        # The speed of the treadmill belt for the planted foot is added to the distance
        # between successive steps e.g. For the left step, the right foot is planted,
        # therefore the effective AP distance traveled by the left foot is speed of the right
        # treadmill belt plus the absolute distance between the right and left footstrikes
        stepcoords = diff(getindex.(steps, :pos)) .+
                        ( isodd(i) ? SVector{2}(0.0, 0.96*steptimes[i]) : SVector{2}(0.0, 1.2*steptimes[i])
                         for i in 1:(length(steps)-1) )
    else
        stepcoords = diff(getindex.(steps, :pos)) .+
                        ( SVector{2}(0.0, 1.2*steptimes[i]) for i in 1:(length(steps)-1) )
    end

    # Spatial and temporal step descriptives
    left_stepcoords = stepcoords[1:2:end] # Odd elements
    left_steplavg = mean(abs.(getindex.(left_stepcoords,2)))
    left_steplstd = std(abs.(getindex.(left_stepcoords,2)))

    right_stepcoords = stepcoords[2:2:end] # Even elements
    right_steplavg = mean(abs.(getindex.(right_stepcoords,2)))
    right_steplstd = std(abs.(getindex.(right_stepcoords,2)))

    # Push into results Dict
    results[:left_steplavg] = left_steplavg
    results[:left_steplstd] = left_steplstd
    results[:right_steplavg] = right_steplavg
    results[:right_steplstd] = right_steplstd

    left_avgsteptime = mean(steptimes[1:2:end]) # Odd elements
    right_avgsteptime = mean(steptimes[2:2:end]) # Even elements

    # Asymmetry indices
    results[:stepasym_spatial] = asymmetry(left_steplavg, right_steplavg)
    results[:stepasym_temporal] = asymmetry(left_avgsteptime, right_avgsteptime)

    # Step width isn't confounded by induced asymmetric gait
    results[:stepwavg] = mean(abs.(getindex.(stepcoords, 1)))
    results[:stepwstd] = std(abs.(getindex.(stepcoords, 1)))

    ########################################
    # Arm swing
    ########################################
    rfc = round.(Int, seg.data.events[:RFC][1:numstrides+1] .* fs) .- 1 # Subtract 1 because the θ and ω have the first element trimmed

    # Drop last element (it is a NaN)
    rhip_θ = seg.data.data[:RHip][1:end-1,1]
    lhip_θ = seg.data.data[:LHip][1:end-1,1]
    rsho_θ = seg.data.data[:RShoulder][1:end-1,1]
    lsho_θ = seg.data.data[:LShoulder][1:end-1,1]

    # ROM and swing asymmetry
    rsho_avgrom, rsho_stdrom = meanrange(rsho_θ, rfc)
    lsho_avgrom, lsho_stdrom = meanrange(lsho_θ, rfc)
    results[:swingasym] = asymmetry(lsho_avgrom, rsho_avgrom)

    results[:rsho_avgrom] = rsho_avgrom
    results[:lsho_avgrom] = lsho_avgrom
    results[:rsho_stdrom] = rsho_stdrom
    results[:lsho_stdrom] = lsho_stdrom

    # Continuous phase
    rhip_crp = continuousphase(rhip_θ, rfc)
    lhip_crp = continuousphase(lhip_θ, rfc)
    rsho_crp = continuousphase(rsho_θ, rfc)
    lsho_crp = continuousphase(lsho_θ, rfc)

    # The subtraction of the two angles has the possibility of introducing discontinuties (which
    # will interfere with the interpolation of the time normalizing)
    crp_lsho_rhip = unwrap(rhip_crp; range=2pi) - unwrap(lsho_crp; range=2pi)
    crp_rsho_lhip = unwrap(lhip_crp; range=2pi) - unwrap(rsho_crp; range=2pi)

    norm_crp_lsho_rhip = timenormalize(crp_lsho_rhip, rfc)
    norm_crp_rsho_lhip = timenormalize(crp_rsho_lhip, rfc)

    sdcrp_lsho_rhip = circstd.(eachrow(reshape(norm_crp_lsho_rhip, 100, numstrides)))
    msdcrp_lsho_rhip = circmean(sdcrp_lsho_rhip)

    sdcrp_rsho_lhip = circstd.(eachrow(reshape(norm_crp_rsho_lhip, 100, numstrides)))
    msdcrp_rsho_lhip = circmean(sdcrp_rsho_lhip)

    results[:msdcrp_lsho_rhip] = rad2deg(msdcrp_lsho_rhip)
    results[:msdcrp_rsho_lhip] = rad2deg(msdcrp_rsho_lhip)

    ########################################
    ## Lyapunov exponent
    ########################################
    frames = round(Int, seg.data.events[:RFC][1]*fs):round(Int, seg.data.events[:RFC][numstrides+1]*fs)

    data = Matrix{Float64}(undef, length(frames), 2*3)
    three = 1:3
    for (i, t) in enumerate(["TrunkLinVel", "TrunkAngVel"])
        data[:,three .+ (i-1)*3] = seg.data.data[Symbol(t)][frames,:]
    end

    NW = Matrix{Float64}(undef, 12_500, 2*3)
    for i in 1:size(NW,2)
        itp = interpolate(data[:,i], BSpline(Cubic(Line(Interpolations.OnGrid()))))
        NW[:,i] .= itp(range(1, stop=length(frames), length=12_500))
    end

    linvelstd = sum(std(view(NW, :, 1:3); dims=1))
    angvelstd = sum(std(view(NW, :, 4:6); dims=1))

    NW[:,1:3] ./= linvelstd
    NW[:,4:6] ./= angvelstd

    d = Dataset(NW)
    r = reconstruct(d, 2, 25)

    ks = 0:50
    E = numericallyapunov(r, ks; w=100, ntype=FixedMassNeighborhood(1), distance=Euclidean())

    λ = ([ ones(length(ks)) ks./100 ] \ E)[2]
    results[:lambdaS] = λ

    return AnalyzedSegment(seg, results)
end

end


text:
using LinearAlgebra

@testset "inner products" begin
    q = [ x->1/sqrt(2), x->x*sqrt(3/2), x->(3x^2-1)*sqrt(5/8), x->(5x^3-3x)*sqrt(7/8) ]
    p = LegendreInnerProduct(length(q))
    x = evalpoints(p)
    B = Basis(q)
    Q = [q(x) for x in x, q in q]

    for i = 1:length(q)
        @test @inferred(norm(q[i]    ,p)) ≈ 1
        @test @inferred(norm(q[i].(x),p)) ≈ 1
    end
    for i = 1:length(q), j = 1:length(q)
        @test abs(@inferred(dot(q[i]  ,q[j]  ,p)) - I[i,j]) <= sqrt(eps())
        @test abs(@inferred(dot(Q[:,i],q[j]  ,p)) - I[i,j]) <= sqrt(eps())
        @test abs(@inferred(dot(q[i]  ,Q[:,j],p)) - I[i,j]) <= sqrt(eps())
        @test abs(@inferred(dot(Q[:,i],Q[:,j],p)) - I[i,j]) <= sqrt(eps())
    end
    for i = 1:length(q)
        e = zeros(length(q)); e[i] = 1
        @test dot(q[i],q,p) ≈ e'
        @test dot(q[i],Q,p) ≈ e'
        @test dot(q[i],B,p) ≈ e'
        @test dot(Q[:,i],q,p) ≈ e'
        @test dot(Q[:,i],Q,p) ≈ e'
        @test dot(Q[:,i],B,p) ≈ e'
        @test dot(q,q[i],p) ≈ e
        @test dot(Q,q[i],p) ≈ e
        @test dot(B,q[i],p) ≈ e
        @test dot(q,Q[:,i],p) ≈ e
        @test dot(Q,Q[:,i],p) ≈ e
        @test dot(B,Q[:,i],p) ≈ e
    end
    @test dot(q,q,p) ≈ I
    @test dot(q,Q,p) ≈ I
    @test dot(q,B,p) ≈ I
    @test dot(Q,q,p) ≈ I
    @test dot(Q,Q,p) ≈ I
    @test dot(Q,B,p) ≈ I
    @test dot(B,q,p) ≈ I
    @test dot(B,Q,p) ≈ I
    @test dot(B,B,p) ≈ I
end


text:
# a normal HDF5 dataset is not a subtype of AbstractArray, but then again HDF5.jl does not perform
# ahead-of-time bounds checking. checkbounds() and CartesianIndices are only defined for AbstractArrays
# and I don't want to also reimplement those
struct SparseDataset{T} <: AbstractArray{T, 2}
    group::HDF5.Group
    csr::Bool

    function SparseDataset(group::HDF5.Group)
        csr = read_attribute(group, "encoding-type")[1:3] == "csr"
        return new{eltype(group["data"])}(group, csr)
    end
end

function Base.getproperty(dset::SparseDataset, s::Symbol)
    if s === :id
        return getfield(dset, :group).id
    elseif s === :file
        return getfield(dset, :group).file
    elseif s === :xfer
        return Tuple(d.xfer for d in getfield(dset, :group))
    else
        return getfield(dset, s)
    end
end

HDF5.filename(dset::SparseDataset) = HDF5.filename(dset.group)
HDF5.copy_object(
    src_obj::SparseDataset,
    dst_parent::Union{HDF5.File, HDF5.Group},
    dst_path::AbstractString,
) = copy_object(src_obj.group, dst_parent, dst_path)
HDF5.isvalid(dset::SparseDataset) = isvalid(dset.group)

getcolptr(dset::SparseDataset) = dset.group["indptr"]
rowvals(dset::SparseDataset) = dset.group["indices"]
nonzeros(dset::SparseDataset) = dset.group["data"]

Base.isassigned(dset::HDF5.Dataset, i) = i <= length(dset)

Base.ndims(dset::SparseDataset) = length(attributes(dset.group)["shape"])
Base.size(dset::SparseDataset) = Tuple(read_attribute(dset.group, "shape"))
Base.size(dset::SparseDataset, d::Integer) = size(dset)[d]
Base.length(dset::SparseDataset) = prod(size(dset))

rawsize(dset::SparseDataset) = dset.csr ? reverse(size(dset)) : size(dset)

Base.lastindex(dset::SparseDataset) = length(dset)
Base.lastindex(dset::SparseDataset, d::Int) = size(dset, d)
HDF5.datatype(dset::SparseDataset) = datatype(dset.group["data"])
Base.read(dset::SparseDataset) = read_matrix(dset.group)

Base.getindex(dset::SparseDataset, i::Integer) = getindex(dset, CartesianIndices(dset)[i])
Base.getindex(dset::SparseDataset, I::Tuple{Integer, Integer}) = getindex(dset, I[1], I[2])
function Base.getindex(dset::SparseDataset, i::Integer, j::Integer)
    @boundscheck checkbounds(dset, i, j)
    if dset.csr
        i, j = j, i
    end
    colptr = getcolptr(dset)
    rowval = rowvals(dset)
    c1, c2 = colptr[j] + 1, colptr[j + 1]
    rowidx = findfirst(x -> x == i - 1, rowval[c1:c2])
    return rowidx === nothing ? eltype(dset)(0) : nonzeros(dset)[c1 + rowidx - 1]
end
Base.getindex(dset::SparseDataset, ::Colon, ::Colon) = read(dset)
function Base.getindex(dset::SparseDataset, I::AbstractUnitRange, j::Integer)
    @boundscheck checkbounds(dset, I, j)
    return dset.csr ? _getindex(dset, j, I) : _getindex(dset, I, j)
end
function Base.getindex(dset::SparseDataset, i::Integer, J::AbstractUnitRange)
    @boundscheck checkbounds(dset, i, J)
    return dset.csr ? _getindex(dset, J, i) : _getindex(dset, i, J)
end

function Base.getindex(dset::SparseDataset, I::AbstractUnitRange, J::AbstractUnitRange)
    @boundscheck checkbounds(dset, I, J)
    sz = size(dset)
    if first(I) == 1 && last(I) == sz[1] && first(J) == 1 && last(J) == sz[2]
        return dset[:, :]
    end
    if dset.csr
        I, J = J, I
    end
    colptr = getcolptr(dset)
    rows = rowvals(dset)
    data = nonzeros(dset)

    newcols = Vector{eltype(colptr)}(undef, length(J) + 1)
    newcols[1] = 1
    newrows = Vector{eltype(rows)}()
    newdata = Vector{eltype(dset)}()
    for (nc, c) in enumerate(J)
        c1, c2 = colptr[c] + 1, colptr[c + 1]
        currrows = rows[c1:c2] .+ convert(eltype(rows), 1)
        rowidx = findall(x -> x ∈ I, currrows)
        newcols[nc + 1] = newcols[nc] + length(rowidx)

        if length(rowidx) > 0
            currdata = data[c1:c2][rowidx]
            currrows = currrows[rowidx]
            sort!(rowidx, currdata, currrows)
            append!(newrows, currrows .- convert(eltype(newrows), first(I)) .+ convert(eltype(newrows), 1))
            append!(newdata, currdata)
        end
    end
    mat = SparseMatrixCSC(length(I), length(J), newcols, newrows, newdata)
    return dset.csr ? mat' : mat
end

function Base.getindex(dset::SparseDataset, i::Integer, ::Colon)
    @boundscheck checkbounds(dset, i, :)
    return dset.csr ? _getindex(dset, :, i) : _getindex(dset, i, :)
end

function Base.getindex(dset::SparseDataset, ::Colon, i::Integer)
    @boundscheck checkbounds(dset, :, i)
    return dset.csr ? _getindex(dset, i, :) : _getindex(dset, :, i)
end

Base.getindex(dset::SparseDataset, I::AbstractUnitRange, ::Colon) = dset[I, 1:size(dset, 2)]

Base.getindex(dset::SparseDataset, ::Colon, J::AbstractUnitRange) = dset[1:size(dset, 1), J]

function _getindex(dset, i::Integer, J::AbstractUnitRange)
    colptr = getcolptr(dset)
    c1, c2 = colptr[first(J)] + 1, colptr[last(J) + 1]
    rows = rowvals(dset)[c1:c2] .+ 1
    rowidx = findall(x -> x == i, rows)

    if length(rowidx) == 0
        return SparseVector(length(J), Vector{eltype(rowidx)}(), Vector{eltype(dset)}())
    end

    data = nonzeros(dset)[c1:c2][rowidx]

    cols = Vector{eltype(rowidx)}(undef, length(data))
    i = 1
    for j in J
        if colptr[j] < rowidx[i] + c1 - 1 <= colptr[j + 1]
            cols[i] = j - first(J) + 1
            if i == length(rowidx)
                break
            end
            i += 1
        end
    end

    sort!(rowidx, data)
    return SparseVector(length(J), cols, data)
end

function _getindex(dset, I::AbstractUnitRange, j::Integer)
    colptr = getcolptr(dset)
    c1, c2 = colptr[j] + 1, colptr[j + 1]
    rows = rowvals(dset)[c1:c2] .+ 1
    rowidx = findall(x -> x ∈ I, rows)
    data = nonzeros(dset)[c1:c2][rowidx]

    sort!(rowidx, data)
    return SparseVector(length(I), rows[rowidx] .- first(I) .+ 1, data)
end

function _getindex(dset, i::Integer, ::Colon)
    rows = read(rowvals(dset))
    nz = nonzeros(dset)
    rowidx = findall(x -> x + 1 == i, rows)
    data = [nz[i] for i in rowidx]

    cols = Vector{eltype(rowidx)}(undef, length(data))
    colptr = getcolptr(dset)
    i = 1
    for (cstart, cend) in zip(1:(length(colptr) - 1), 2:length(colptr))
        if colptr[cstart] < rowidx[i] <= colptr[cend]
            cols[i] = cstart
            if i == length(rowidx)
                break
            end
            i += 1
        end
    end

    return SparseVector(rawsize(dset)[2], cols, data)
end

function _getindex(dset, ::Colon, j::Integer)
    colptr = getcolptr(dset)
    rows = read(rowvals(dset))
    c1, c2 = colptr[j] + 1, colptr[j + 1]
    rowidx = rows[c1:c2] .+ 1
    data = nonzeros(dset)[c1:c2]

    sort!(rowidx, data)
    return SparseVector(rawsize(dset)[1], rowidx, data)
end

function Base.setindex!(dset::SparseDataset, x::Number, i::Integer, j::Integer)
    @boundscheck checkbounds(dset, i, j)
    if dset.csr
        i, j = j, i
    end
    cols = getcolptr(dset)
    rows = rowvals(dset)

    c1, c2 = cols[j] + 1, cols[j + 1]
    rowidx = findfirst(x -> x + 1 == i, rows[c1:c2])
    if rowidx === nothing && x != 0
        throw(KeyError("changing the sparsity structure of a SparseDataset is not supported"))
    elseif x != 0
        nonzeros(dset)[c1 + rowidx - 1] = x
    end
end
function Base.setindex!(
    dset::SparseDataset{<:Number},
    x::AbstractArray{<:Number, 2},
    I::AbstractUnitRange,
    J::AbstractUnitRange,
)
    @boundscheck checkbounds(dset, I, J)
    length(x) == length(I) * length(J) || throw(
        DimensionMismatch(
            "tried to assign $(length(x)) elements to destination of size $(length(I) * length(J))",
        ),
    )
    if dset.csr
        I, J = J, I
        x = x'
    end
    linxidx = LinearIndices(x)
    cols = getcolptr(dset)
    rows = rowvals(dset)

    xidx = Int[]
    dsetidx = Int[]
    for (ic, c) in enumerate(J)
        c1, c2 = cols[c] + 1, cols[c + 1]
        crows = rows[c1:c2] .+ 1
        rowidx = findall(x -> x ∈ I, crows)
        xvals = x[I[I .∉ ((@view crows[rowidx]),)] .- first(I) .+ 1, ic]
        if length(rowidx) != length(I) && any(xvals .!= 0)
            throw(KeyError("changing the sparsity structure of a SparseDataset is not supported"))
        end
        append!(xidx, linxidx[crows[rowidx] .- first(I) .+ 1, ic])
        append!(dsetidx, c1 - 1 .+ rowidx)
    end
    # HDF5 doesn't support assignment using Arrays as indices, so we have to loop
    nz = nonzeros(dset)
    for (di, dx) in zip(dsetidx, @view x[xidx])
        nz[di] = dx
    end
end

Base.eachindex(dset::SparseDataset) = CartesianIndices(size(dset))

function Base.show(io::IO, dset::SparseDataset)
    if isvalid(dset)
        print(io, "Sparse HDF5 dataset: ", HDF5.name(dset.group), " (file: ", dset.file.filename, " xfer_mode: ", Tuple(x.id for x in dset.xfer), ")")
    else
        print(io, "Sparse HDF5 datset: (invalid)")
    end
end

function Base.show(io::IO, ::MIME"text/plain", dset::SparseDataset)
    if get(io, :compact, false)::Bool
        show(io, dset)
    else
        print(io, HDF5._tree_icon(HDF5.Dataset), " ", dset)
    end
end


text:
using Polymer
using Documenter

makedocs(;
    modules=[Polymer],
    authors="Yi-Xin Liu <lyx@fudan.edu.cn> and contributors",
    repo="https://github.com/liuyxpp/Polymer.jl/blob/{commit}{path}#L{line}",
    sitename="Polymer.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://liuyxpp.github.io/Polymer.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/liuyxpp/Polymer.jl",
)


text:
@testset "clustering" begin
    @testset "partitional clustering" begin
        r = [10, 20]
        c = [2, 3, 5]
        C = [0 0 0; 0 0 0; 0 0 0]
        W = [0 0 0; 0 0 0; 0 0 0]
        Y = [1 0 0.5; 0 1 0.5]
        p = (μ=[0 1.0; 1.0 0],)
        pc = PartitionalClustering(r, c, C, W, Y, p)
        pc2 = PartitionalClustering(r, c, C, W, Y, p)
        pc3 = PartitionalClustering(r, c, C, W, Y, p)

        @testset "constructors" begin
            @test isa(pc, PartitionalClustering)
            @test (pc.r == r && pc.c == c && pc.C == C && pc.W == W && pc.Y == Y
                    && pc.p == p)
        end

        @testset "interface constructors" begin
            pc4 = PartitionalClustering(KmeansResult(p.μ, [1, 2, 1],
                    [0.5, 0.25, 0.125], [2, 1], [2, 1], 0.875, 100, true))
            @test isa(pc4, PartitionalClustering)
            @test (pc4.r == Int[] && pc4.c == Int[]
                    && pc4.C == Matrix{Int}(undef, 0, 0)
                    && pc4.W == Matrix{Float64}(undef, 0, 0)
                    && pc4.Y == [1 0 1; 0 1 0]
                    && pc4.p == (centers=p.μ, costs=[0.5, 0.25, 0.125],
                            counts=[2, 1], wcounts=[2, 1], totalcost=0.875,
                            iterations=100, converged=true))
            pc4 = PartitionalClustering(KmedoidsResult([11, 23], [1, 2, 1],
                    [0.5, 0.25, 0.125], [2, 1], 0.875, 100, true))
            @test isa(pc4, PartitionalClustering)
            @test (pc4.r == Int[] && pc4.c == Int[]
                    && pc4.C == Matrix{Int}(undef, 0, 0)
                    && pc4.W == Matrix{Float64}(undef, 0, 0)
                    && pc4.Y == [1 0 1; 0 1 0]
                    && pc4.p == (medoids=[11, 23], costs=[0.5, 0.25, 0.125],
                            counts=[2, 1], totalcost=0.875, iterations=100,
                            converged=true))
            pc4 = PartitionalClustering(FuzzyCMeansResult(p.μ,
                    [1 0; 0.5 0.5; 1 0], 100, true))
            @test isa(pc4, PartitionalClustering)
            @test (pc4.r == Int[] && pc4.c == Int[]
                    && pc4.C == Matrix{Int}(undef, 0, 0)
                    && pc4.W == Matrix{Float64}(undef, 0, 0)
                    && pc4.Y == [1 0.5 1; 0 0.5 0]
                    && pc4.p == (centers=p.μ, iterations=100, converged=true))
        end

        @testset "equality operator" begin
            @test pc == pc
            @test pc == pc2 && pc2 == pc
            @test pc == pc2 && pc2 == pc3 && pc == pc3
        end

        @testset "hash" begin
            @test hash(pc) == hash(pc)
            @test pc == pc2 && hash(pc) == hash(pc2)
        end

        @testset "accessors" begin
            @test axes(pc) == (r, c)
            @test features(pc) == axes(pc, 1) == r
            @test instances(pc) == axes(pc, 2) == c
            @test constraints(pc) == C
            @test weights(pc) == W
            @test parameters(pc) == p
            @test assignments(pc) == Y
        end
    end

    @testset "hierarchical clustering" begin
        r = [10, 20]
        c = [2, 3, 5]
        C = rand([-1, 0, 1], 3, 3, 3)
        W = rand(3, 3, 3) .+ 1
        p = NamedTuple()
        hc = HierarchicalClustering(r, c, C, W, p)

        @testset "constructors" begin
            @test isa(hc, HierarchicalClustering)
            @test hc.r == r && hc.c == c && hc.C == C && hc.W == W && hc.p == p
        end

        @testset "accessors" begin
            @test axes(hc) == (r, c)
            @test features(hc) == axes(hc, 1) == r
            @test instances(hc) == axes(hc, 2) == c
            @test constraints(hc) == C
            @test weights(hc) == W
            @test parameters(hc) == p
        end
    end
end