text:
abstract Constraint
abstract ConvexConstraint <: Constraint
immutable BoxConstraint <: ConvexConstraint
    lower
    upper
end

immutable BallConstraint <: ConvexConstraint
    center # center
    radius # radius
    p # l^p norm
end

immutable EqualityConstraint <: Constraint
    c
    Jc!
    cJc!
end

function EqualityConstraint(c, Jc!)
    function cJc!(x::Array, storage::Array)
        Jc!(x, storage)
        return c(x)
    end
    return EqualityConstraint(c, Jc!, cJc!)
end

type AugmentedLagrangian{Tl<:Union{Real, Vector}}
   F::DifferentiableFunction
   C::EqualityConstraint
   lambda::Tl
   mu::Real
   Dc::Matrix
end

typealias AL AugmentedLagrangian
function AugmentedLagrangian(F::DifferentiableFunction,
                 C::EqualityConstraint,
                 x0::AbstractVector;
                 lambda=:auto, mu=10.0)
   # TODO: maybe not ideal to evaluate F, C here, maybe better to set
   # dimensions when C is first called??? But see #282
   C0 = C.c(x0)
   dim_x = length(x0)
   dim_c = length(C0)
   if lambda == :auto
      lambda = - mu * C0
   end
   return AugmentedLagrangian(F, C, lambda, mu, zeros(dim_c, dim_x))
end

max_stars_repo_name: cortner/AugmentedLagrangian

max_issues_repo_name: cortner/AugmentedLagrangian

max_forks_repo_name: cortner/AugmentedLagrangian


text:
#
# ## Pendulum
# # https://github.com/HarvardAgileRoboticsLab/unscented-dynamic-programming/blob/master/pendulum_dynamics.m
# function pendulum_dynamics!(ẋ::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T}) where T
#     m = 1.
#     l = 0.5
#     b = 0.1
#     lc = 0.5
#     I = 0.25
#     g = 9.81
#     ẋ[1] = x[2]
#     ẋ[2] = (u[1] - m*g*lc*sin(x[1]) - b*x[2])/I
# end
#
# n,m = 2,1
# pendulum = Model(pendulum_dynamics!,n,m) # inplace model
#
# # unknown mass
# function pendulum_dynamics_uncertain!(ẋ,x,u,w)
#     m = 1. + w[1]
#     l = 0.5
#     b = 0.1
#     lc = 0.5
#     I = 0.25
#     g = 9.81
#
#     ẋ[1] = x[2]
#     ẋ[2] = u[1]/(m*lc*lc) - g*sin(x[1])/lc - b*x[2]/(m*lc*lc)
#     return nothing
# end
#
# n = 2; m = 1; r = 1
# pendulum_uncertain = UncertainModel(pendulum_dynamics_uncertain!,n,m,r)

@with_kw mutable struct Pendulum{T} <: AbstractModel
    mass::T = 1.
    length::T = 0.5
    b::T = 0.1
    lc::T = 0.5
    I::T = 0.25
    g::T = 9.81
end

Base.size(::Pendulum) = 2,1

function dynamics(p::Pendulum, x, u)
    m = p.mass * p.lc * p.lc
    @SVector [x[2],
              u[1]/m - p.g*sin(x[1])/p.lc - p.b*x[2]/m]
end
@inline Base.position(::Pendulum, x) = @SVector zeros(3)
orientation(::Pendulum, x) = expm((pi-x[1])*@SVector [1,0,0.])

ext: jl


text:
"""
random vector on unit sphere, should be uniformly distributed?
"""
function rand_sphere(T = Float64)
   R = randn(SVector{3, T})
   return R / norm(R)
end

"""
random rotation matrix;
WARN: never checked what the distribution is
"""
rand_rot() = (K = (@SMatrix rand(3,3)) .- 0.5; exp(K - K'))

"""
random reflection, represented as 1, -1 Integer
"""
rand_refl() = rand([-1,1])

"""
random isometry i.e. element of O(3)
"""
rand_O3() = rand_refl() * rand_rot()

# the following two functions are a little hack to make sure
# that the basis spec is read in the same symbol-order as it is written
# (since dicts don't have a specified ordering...)

function _write_dict_1pspec(spec::Vector{NamedTuple{SYMS, NTuple{NSYM, Int}}}) where {SYMS, NSYM}
   inds = Vector{Int}[]
   for b in spec
      vals = [getproperty(b, sym) for sym in SYMS]
      push!(inds, vals)
   end
   return Dict("SYMS" => [ string.(SYMS)... ], "inds" => inds)
end

function _read_dict_1pspec(D::Dict)
   NTPROTO = namedtuple(D["SYMS"]...)
   return [ NTPROTO(binds) for binds in D["inds"] ]
end


text:
# Testing discrete univariate distributions

using Distributions
import JSON
using Base.Test
using Compat

function verify_and_test_drive(jsonfile, selected, n_tsamples::Int)
    R = JSON.parsefile(jsonfile)
    for (ex, dct) in R
        dsym = Symbol(dct["dtype"])
        dname = string(dsym)

        # test whether it is included in the selected list
        # or all are selected (when selected is empty)
        if !isempty(selected) && !(dname in selected)
            continue
        end

        # perform testing
        println("    testing $(ex)")
        dtype = eval(dsym)
        d = eval(parse(ex))
        if dtype == TruncatedNormal
            @test isa(d, Truncated{Normal{Float64}})
        else
            @assert isa(dtype, Type) && dtype <: UnivariateDistribution
            @test isa(d, dtype)
        end

        # verification and testing
        verify_and_test(d, dct, n_tsamples)
    end
end

_parse_x(d::DiscreteUnivariateDistribution, x) = round(Int, x)
_parse_x(d::ContinuousUnivariateDistribution, x) = Float64(x)

_json_value(x::Number) = x

_json_value(x::AbstractString) =
    x == "inf" ? Inf :
    x == "-inf" ? -Inf :
    x == "nan" ? NaN :
    error("Invalid numerical value: $x")

function verify_and_test(d::UnivariateDistribution, dct::Dict, n_tsamples::Int)
    # verify parameters
    pdct = dct["params"]
    for (fname, val) in pdct
        f = eval(Symbol(fname))
        @assert isa(f, Function)
        Base.Test.test_approx_eq(f(d), val, "$fname(d)", "val")
    end

    # test various constructors for promotion, all-Integer args, etc.
    D = eval(Symbol(dct["dtype"]))  # distribution name
    pars = params(d)

    # promotion constructor:
    float_pars = map(x -> isa(x, AbstractFloat), pars)
    if length(pars) > 1 && sum(float_pars) > 1
        mixed_pars = Any[pars...]
        first_float = findfirst(float_pars)
        mixed_pars[first_float] = Float32(mixed_pars[first_float])

        @test typeof(D(mixed_pars...)) == typeof(d)
    end

    # promote integer arguments to floats, where applicable
    if sum(float_pars) >= 1 && !any(map(isinf, pars)) && !isa(d, Geometric)
        int_pars = map(x -> ceil(Int, x), pars)
        @test typeof(D(int_pars...)) == typeof(d)
    end

    # verify stats
    @test_approx_eq minimum(d) _json_value(dct["minimum"])
    @test_approx_eq maximum(d) _json_value(dct["maximum"])
    @test_approx_eq_eps mean(d) _json_value(dct["mean"]) 1.0e-8
    if !isa(d, VonMises)
        @test_approx_eq_eps var(d) _json_value(dct["var"]) 1.0e-8
    end
    if !isa(d, Skellam)
        @test_approx_eq_eps median(d) _json_value(dct["median"]) 1.0
    end

    if applicable(entropy, d) && !isa(d, VonMises)  # SciPy VonMises entropy is wrong
        @test_approx_eq_eps entropy(d) dct["entropy"] 1.0e-7
    end

    # test conversions if distribution is parametric
    if !isempty(typeof(d).parameters) && !isa(d, Truncated)
        D = typeof(d).name.primary
        W = Float32
        @test typeof(convert(D{W}, d)) == D{W}
        @test typeof(convert(D{W}, params(d)...)) == D{W}
    end

    # verify quantiles
    if !isa(d, Union{Skellam, VonMises})
        @test_approx_eq_eps quantile(d, 0.10) dct["q10"] 1.0e-8
        @test_approx_eq_eps quantile(d, 0.25) dct["q25"] 1.0e-8
        @test_approx_eq_eps quantile(d, 0.50) dct["q50"] 1.0e-8
        @test_approx_eq_eps quantile(d, 0.75) dct["q75"] 1.0e-8
        @test_approx_eq_eps quantile(d, 0.90) dct["q90"] 1.0e-8
    end

    # verify logpdf and cdf at certain points
    pts = dct["points"]
    for pt in pts
        x = _parse_x(d, pt["x"])
        lp = Float64(pt["logpdf"])
        cf = Float64(pt["cdf"])
        Base.Test.test_approx_eq(logpdf(d, x), lp, "logpdf(d, $x)", "lp")
        if !isa(d, Skellam)
            Base.Test.test_approx_eq(cdf(d, x), cf, "cdf(d, $x)", "cf")
        else
            @test_throws MethodError cdf(d, x)
        end
    end

    try
        m = mgf(d,0.0)
        @test m == 1.0
    catch e
        isa(e, MethodError) || throw(e)
    end
    try
        c = cf(d,0.0)
        @test c == 1.0
        # test some extra values: should all be well-defined
        for t in (0.1,-0.1,1.0,-1.0)
            @test !isnan(cf(d,t))
        end
    catch e
        isa(e, MethodError) || throw(e)
    end

    # generic testing
    if isa(d, Cosine)
        n_tsamples = floor(Int, n_tsamples / 10)
    end
    if !isa(d, Union{Skellam, VonMises})
        test_distr(d, n_tsamples)
    end
end

## main

for c in ["discrete",
          "continuous",
          "discrete_hand_coded"]

    title = string(uppercase(c[1]), c[2:end])
    println("    [$title]")
    println("    ------------")
    jsonfile = joinpath(dirname(@__FILE__), "$(c)_test.json")
    verify_and_test_drive(jsonfile, ARGS, 10^6)
    println()
end


text:
module TestDDMortality
    using BioEnergeticFoodWebs
    using Test

    A = [0 1 0 ; 0 0 1 ; 0 0 0]
    ddm = (x -> x .* 0.2)
    p = model_parameters(A, dc = ddm)

    b = [0.5, 0.6, 0.8]
    PI_death = BioEnergeticFoodWebs.density_dependent_mortality(p, b)
    @test all(PI_death .== ddm(b) .* Int.(.!p[:is_producer]))

    dbdt = BioEnergeticFoodWebs.dBdt(zeros(3), b, p, 1)
    # Consumption
    gain, loss = BioEnergeticFoodWebs.consumption(p, b)
    # Growth
    growth, G = BioEnergeticFoodWebs.get_growth(p, b)
    # Balance
    balance = zeros(eltype(b), length(b))
    for i in eachindex(balance)
      balance[i] = growth[i] + gain[i] - loss[i] - PI_death[i]
    end

    @test all(dbdt .== balance)

    ddm_prod = (x -> x .^ 2 .* 0.1)
    p = model_parameters(A, dc = ddm, dp = ddm_prod)
    PI_death = BioEnergeticFoodWebs.density_dependent_mortality(p, b)
    mc = ddm(b).* Int.(.!p[:is_producer])
    mp = ddm_prod(b).* Int.(p[:is_producer])
    @test all(PI_death .== mc .+ mp)

    dbdt = BioEnergeticFoodWebs.dBdt(zeros(3), b, p, 1)
    # Consumption
    gain, loss = BioEnergeticFoodWebs.consumption(p, b)
    # Growth
    growth, G = BioEnergeticFoodWebs.get_growth(p, b)
    # Balance
    balance = zeros(eltype(b), length(b))
    for i in eachindex(balance)
      balance[i] = growth[i] + gain[i] - loss[i] - PI_death[i]
    end

    @test all(dbdt .== balance)

end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
using Documenter, TidalFluxQuantities

makedocs()

deploydocs(deps=Deps.pip("mkdocs","python-markdown-math"),
           repo="github.com/wkearn/TidalFluxQuantities.jl.git",
           julia = "release"
)


text:
@testset "Matrices" begin
  # tests with Aitchson's definitions
  d = 10
  D = d+1
  @test J(d) == ones(d, d)
  @test F(d) == [I(d) -ones(d)]
  @test G(D) ≈ I(D) - J(D)/ D
  @test H(d) == I(d) + J(d)

  # tests of multiplications
  for M in [J G H]
    @test M * [1:d...] ≈ M(d) * [1:d...]
    @test [1:d...]' * M ≈ [1:d...]' * M(d)
  end

  # tests of F (d x D)
  @test F * [1:D...] == F(d) * [1:D...]
  @test [1:d...]' * F == [1:d...]' * F(d)
end


text:
struct FPFState{T, TF} <: AbstractFilterState
    ensemble::UnweightedParticleEnsemble{T}
    eq::PoissonEquation{T, TF}
    function FPFState(ens::UnweightedParticleEnsemble{T}, eq::PoissonEquation{T, TF}) where {T, TF}
        update!(eq, ens)
        new{T, TF}(ens, eq)
    end
end

struct FPF{TS<:HiddenStateModel, TO<:ObservationModel, T<:GainEstimationMethod} <: AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime}
    state_model::TS
    obs_model::TO
    gain_method::T
    N::Int
    function FPF(st_mod::HiddenStateModel{S1, ContinuousTime}, ob_mod::ObservationModel{S1, S2, ContinuousTime}, gain_method::GainEstimationMethod, N::Int) where {S1, S2}
        if ob_mod isa DiffusionObservationModel || ob_mod isa LinearDiffusionObservationModel
            return new{typeof(st_mod), typeof(ob_mod), typeof(gain_method)}(st_mod, ob_mod, gain_method, N)
        else
            error("Error: the FPF algorithm is not defined for non-diffusion-type observations. Try ppFPF for counting process observations.")
        end
    end
end

#####################
###### METHODS ######
#####################

initial_condition(fpf::FPF) = fpf.state_model.init
no_of_particles(fpf::FPF) = fpf.N
no_of_particles(st::FPFState) = no_of_particles(st.ensemble)
state_model(fpf::FPF) = fpf.state_model
obs_model(fpf::FPF) = fpf.obs_model
gain_estimation_method(fpf::FPF) = fpf.gain_method

mean(st::FPFState) = mean(st.ensemble)
cov(st::FPFState)  = cov(st.ensemble)
var(st::FPFState)  = var(st.ensemble)

function Base.show(io::IO, ::MIME"text/plain", fpf::FPF)
    print(io, "Feedback particle filter algorithm
    hidden state model:                     ", state_model(fpf),"
    observation model:                      ", obs_model(fpf),"
    ensemble size:                          ", no_of_particles(fpf),"
    gain estimation method:                 ", gain_estimation_method(fpf))
end

function Base.show(io::IO, fpf::FPF)
    print(io, "FPF with ", no_of_particles(fpf)," particles and ", gain_estimation_method(fpf))
end

######################
### MAIN ALGORITHM ###
######################

function initialize(fpf::FPF)
    return FPFState(fpf)
end

function update!(filter_state::FPFState, filt_algo::FPF, dY, dt)
    propagate!(filter_state, filt_algo, dt)
    assimilate!(dY, filter_state, filt_algo, dt)
end

function propagate!(filter_state::FPFState, filt_algo::FPF, dt)
    propagate!(filter_state.ensemble, filt_algo.state_model, dt)
end

function assimilate!(dY, filter_state::FPFState, filt_algo::FPF, dt)
    ensemble  = filter_state.ensemble
    eq        = filter_state.eq
    method    = filt_algo.gain_method

    error     = dY .- dt * eq.H/2 .- dt * eq.mean_H/2
    solve!(eq, method)

    heun!(eq, ensemble, error, method) # stochastic Heun method: intermediate step

    applygain!(ensemble, eq, error)
    update!(eq, ensemble)
    return ensemble.positions
end

########################
### HELPER FUNCTIONS ###
########################

function heun!(eq::PoissonEquation, ensemble::UnweightedParticleEnsemble, error::AbstractMatrix, method::GainEstimationMethod)
    ensemble2 = deepcopy(ensemble)
    applygain!(ensemble2, eq, error)
    eq2       = deepcopy(eq)
    update!(eq2, ensemble2)
    solve!(eq2, method)
    eq.gain  .= (eq.gain + eq2.gain) / 2
end

function gainxerror(gain::Array{T, 3}, error::Array{T, 2}) where T
    out = zeros(T, size(gain, 1), size(gain, 2))
    add_gainxerror!(out, gain, error)
    return out
end

function add_gainxerror!(out::Array{T, 2}, gain::Array{T, 3}, error::Array{T, 2}) where T
    size(gain, 2) == size(error, 2) == size(out, 2) ? N = size(gain, 2) : throw(DimensionMismatch("The provided gain, error, and output array are for different numbers of particles."))
    size(gain, 3) == size(error, 1)                 ? m = size(gain, 3) : throw(DimensionMismatch("The provided gain and error are for different numbers of observed variables."))
    size(gain, 1) == size(out, 1)                   ? n = size(gain, 1) : throw(DimensionMismatch("The provided gain and output array are for different numbers of hidden variables."))
    @inbounds for k in 1:N, i in 1:n, j in 1:m
            out[i, k] += gain[i, k, j] * error[j, k]
    end
    return out
end

function applygain!(ens::UnweightedParticleEnsemble, eq::GainEquation, error::AbstractMatrix)
    add_gainxerror!(ens.positions, eq.gain, error)
end

################################
### CONVENIENCE CONSTRUCTORS ###
################################

function FPFState(fpf::FPF)
    N   = no_of_particles(fpf)
    ens = UnweightedParticleEnsemble(hcat([initialize(fpf.state_model) for i in 1:N]...))
    eq  = PoissonEquation(observation_function(fpf.obs_model), ens)
    return FPFState(ens, eq)
end

function FPFState(filt_prob::AbstractFilteringProblem, N)
    return FPFState(state_model(filt_prob), obs_model(filt_prob), N)
end

function FPFState(state_model::HiddenStateModel, obs_model::ObservationModel, N)
    ens = UnweightedParticleEnsemble(hcat([initialize(state_model) for i in 1:N]...))
    eq  = PoissonEquation(observation_function(obs_model), ens)
    return FPFState(ens, eq)
end

function FPF(filt_prob::AbstractFilteringProblem, method::GainEstimationMethod, N)
    if state_model(filt_prob) isa DiffusionStateModel
        st_mod = state_model(filt_prob)
    elseif state_model(filt_prob) isa LinearDiffusionStateModel
        st_mod = state_model(filt_prob)     # convert linear diffusion model into generic one
    else
        error("Error: cannot build feedback particle filter for given state model. Must be a diffusion.")
    end

    if obs_model(filt_prob) isa DiffusionObservationModel
        ob_mod = obs_model(filt_prob)
    elseif obs_model(filt_prob) isa LinearDiffusionObservationModel
        ob_mod = obs_model(filt_prob)       # convert linear diffusion model into generic one
    else
        error("Error: cannot build feedback particle filter for given observation model. Must be a diffusion.")
    end

    return FPF(st_mod, ob_mod, method, N)
end


text:
prob_clust(P, xhat, cX, yhat, cY) = sum(P[cX.==xhat,cY.==yhat])

function prob_clust(P, Xhat::UnitRange, cX, Yhat::UnitRange, cY)
	output = Array(Float64, (length(Xhat), length(Yhat)))
    for xhat in Xhat
        for yhat in Yhat
            output[xhat, yhat] = prob_clust(P, xhat,cX, yhat,cY)
    	end
	end
	return output
end

prob_x_given_xhat(P, x, xhat, cX) = sum(P[x,:]) / sum(P[cX.==xhat,:])
prob_y_given_yhat(P ,y, yhat, cY) = sum(P[:,y]) / sum(P[:,cY.==yhat])

prob_Y_given_x(P, x) = P[x,:] / sum(P[x,:])
prob_X_given_y(P, y) = P[:,y] / sum(P[:,y])

calc_q(p, x,cX, y,cY) = prob_clust(p, cX[x],cX, cY[y], cY) * prob_x_given_xhat(p, x, cX[x],cX) * prob_y_given_yhat(p, y, cY[y],cY)

function calc_q(p, X::UnitRange, cX, Y::UnitRange, cY)
	output = Array(Float64, (length(X), length(Y)))
	for x in X
        for y in Y
            output[x,y] = calc_q(p, x,cX, y,cY)
    	end
	end
	return output
end

prob_Y_given_xhat(P, xhat,cX) = sum(P[cX.==xhat,:] / sum(P[cX.==xhat,:]), 1)
prob_X_given_yhat(P, yhat,cY) = sum(P[:,cY.==yhat] / sum(P[:,cY.==yhat]), 2)

function next_cx(p,q, x::Int,cX::Array, k)
    q_dist_xhat = Array(Float64, k)
    p_dist_x = prob_Y_given_x(p,x)

    for xhat in 1:k
        q_dist_xhat[xhat] = kl_divergence( vec(p_dist_x), vec(prob_Y_given_xhat(q, xhat,cX)) )
    end
    return indmin(q_dist_xhat)
end

function next_cX(p,q, cX::Array, k)
    output = Array(Int, size(cX))
    for x in 1:size(cX,2)
        output[x] = next_cx(p,q, x,cX, k)
    end
    return output
end

function next_cy(p,q, y::Int,cY::Array{Int,2}, l)
    q_dist_yhat = Array(Float64,l)
    p_dist_y = prob_X_given_y(p,y)
    for yhat in 1:l
        q_dist_yhat[yhat] = kl_divergence( vec(p_dist_y), vec(prob_X_given_yhat(q, yhat, cY)) )
    end
    return indmin(q_dist_yhat)
end

function next_cY(p,q, cY::Array{Int,2}, l)
    output = Array(Int, size(cY))
    for y in 1:size(cY,2)
        output[y] = next_cy(p,q, y, cY, l)
    end
    return output
end

function printState()
    println("cX:",cX, " cY:",cY)
    println(round( q, 3))
    println(prob_clust(p, 1:k,cX, 1:l,cY))
end

max_stars_repo_name: slongwell/ITCC

max_issues_repo_name: slongwell/ITCC

max_forks_repo_name: slongwell/ITCC


text:
## Script for plots of fixed parameters

using Dates, Statistics
using StatsPlots, JLD2, Plots.PlotMeasures, CSV, DataFrames
using Parameters, Distributions, LinearAlgebra
using RecursiveArrayTools, MCMCChains, Interpolations, Suppressor
import KenyaCoVaccines

## Age dependent rates fixed in model
# 0-19, 20-49, 50-59, 60-69, 70-79 and 80+ year olds

model = load("modelfits/Nairobi_model.jld2")["model"]
age_strs = ["0-19", "20-49", "50-59", "60-69", "70-79", "80+"]
rel_sus_plt = bar(model.p_sus,
    lab = "",
    ylabel = "Rel. susceptibility",
    xlabel = "Age group",
    xticks = (1:6, age_strs),
    title = "Relative susceptibility in transmission model",
    tickfont = 11, guidefont = 13, titlefont = 16)

symptom_plt = bar(model.p_symp,
    lab = "",
    ylabel = "Prob. symptoms",
    xlabel = "Age group",
    xticks = (1:6, age_strs),
    title = "Probability of symptomatic disease in transmission model",
    tickfont = 11, guidefont = 13, titlefont = 12)

# savefig(rel_sus_plt, "generation_of_data_for_econ_analysis/SI_plots/rel_sus_plt.png")
# savefig(symptom_plt, "generation_of_data_for_econ_analysis/SI_plots/symp_plt.png")

## Age dependent contact rates fixed in model
max_mixing_rate = maximum([maximum(model.M_county_ho), maximum(model.M_county_other), maximum(model.M_county_school), maximum(model.M_county_work)])

plt_M_home = heatmap(model.M_county_ho,
    clims = (0, max_mixing_rate),
    yticks = (1:6, age_strs),
    xticks = (1:6, age_strs),
    ylabel = "Age group (from)",
    xlabel = "Age group (to)",
    title = "Scaled contact rates from person to group (Home)")

plt_M_other = heatmap(model.M_county_other,
    clims = (0, max_mixing_rate),
    yticks = (1:6, age_strs),
    xticks = (1:6, age_strs),
    ylabel = "Age group (from)",
    xlabel = "Age group (to)",
    title = "Scaled contact rates from person to group (Other)")

plt_M_school = heatmap(model.M_county_school,
    clims = (0, max_mixing_rate),
    yticks = (1:6, age_strs),
    xticks = (1:6, age_strs),
    ylabel = "Age group (from)",
    xlabel = "Age group (to)",
    title = "Scaled contact rates from person to group (School)")

plt_M_work = heatmap(model.M_county_work,
    clims = (0, max_mixing_rate),
    yticks = (1:6, age_strs),
    xticks = (1:6, age_strs),
    ylabel = "Age group (from)",
    xlabel = "Age group (to)",
    title = "Scaled contact rates from person to group (Work)")

# savefig(plt_M_home, "generation_of_data_for_econ_analysis/SI_plots/home_mixing.png")
# savefig(plt_M_other, "generation_of_data_for_econ_analysis/SI_plots/other_mixing.png")
# savefig(plt_M_school, "generation_of_data_for_econ_analysis/SI_plots/school_mixing.png")
# savefig(plt_M_work, "generation_of_data_for_econ_analysis/SI_plots/work_mixing.png")

## Plots for the infection model

@load("data/fitted_mort_by_county_age_variant.jld2")
@load("data/fitted_rel_risk_ICU.jld2")
@load("data/fitted_rel_risk_hosp.jld2")

# countynames
counties_with_serodata = ["Nairobi", "Mombasa",
    "Nakuru", "Uasin Gishu", "Embu", "Kisumu", "Siaya", "Kisii", "Nyeri",
    "Kilifi", "Kwale"]

semiurbancounties_lowserodata = ["Kajiado", "Kiambu", "Machakos", "Isiolo", "Kirinyaga", "Murang'a",
    "Trans Nzoia", "Nyandarua", "Meru", "Bungoma", "Nyamira", "Kakamega",
    "Laikipia", "Taita Taveta", "Makueni", "Elgeyo-Marakwet",
    "Tharaka-Nithi", "Kericho", "Bomet", "Homa Bay", "Busia", "Migori",
    "Nandi", "Vihiga"]

ruralcounties_lowserodata = ["Garissa", "Kitui",
    "Lamu", "Marsabit", "Narok", "Tana River", "Turkana", "Wajir",
    "Baringo", "Mandera", "Samburu", "West Pokot"]

countynames = sort([counties_with_serodata; semiurbancounties_lowserodata; ruralcounties_lowserodata])

plt_rel_report = bar(mort_scales,
    orientation = :horizontal,
    lab = "",
    title = "Relative reporting rate for clinical outcomes",
    yticks = (1:47, countynames),
    xlabel = "Rel. reporting rate",
    xticks = 0:0.1:1, xlims = (0, 1.0),
    tickfont = 9, guidefont = 13, titlefont = 12, dpi = 250,
    size = (500, 900), left_margin = 20mm, right_margin = 20mm)

# savefig(plt_rel_report, "generation_of_data_for_econ_analysis/SI_plots/rel_reporting.png")

plt_mort_age = bar(mort_age,
    lab = "",
    xticks = (1:6, age_strs),
    ylabel = "Risk per infection",
    title = "Baseline fatality risk per infection",
    xlabel = "Age group",
    yticks = 0.0:0.005:0.04,
    tickfont = 9, guidefont = 13, titlefont = 14, dpi = 250)

# savefig(plt_mort_age, "generation_of_data_for_econ_analysis/SI_plots/mort_age.png")
var_str = ["wild-type", "Alpha/Beta variant", "Delta variant"]
plt_mort_var = bar([1.0; mort_var],
    lab = "",
    xticks = (1:3, var_str),
    ylabel = "Rel. risk due to variant",
    title = "Rel. fatality risk by variant",
    xlabel = "Variant",
    tickfont = 9, guidefont = 13, titlefont = 14, dpi = 250)

# savefig(plt_mort_var, "generation_of_data_for_econ_analysis/SI_plots/mort_var.png")

##
plt_mort_age_outcome = groupedbar(mort_age .* [1.0; H_ICU_to_death_fit; H_hosp]',
    lab = ["Deadly" "Critical" "Severe"], legend = :topleft,
    xticks = (1:6, age_strs),
    ylabel = "Risk per infection",
    title = "Baseline risk per infection",
    xlabel = "Age group",
    # yticks = 0.0:0.005:0.04,
    tickfont = 9, guidefont = 13, titlefont = 14, dpi = 250)
savefig(plt_mort_age_outcome, "generation_of_data_for_econ_analysis/SI_plots/mort_age_outcome.png")

plt_mort_var_outcome = groupedbar(hcat([1.0; mort_var], [1.0; mort_var]),
    lab = ["Deadly" "Critical" "Severe"], legend = :topleft,
    # xticks = (1:6, age_strs),
    ylabel = "Risk per infection",
    title = "Baseline risk per infection",
    xlabel = "Age group",
    # yticks = 0.0:0.005:0.04,
    tickfont = 9, guidefont = 13, titlefont = 14, dpi = 250)


text:
aligned, aligned_lm, reco = align_principal(images, lms_on_surface, 192)

resized, resized_lm, scales = resize_relevant(aligned, aligned_lm, 128)
aligned = nothing

images1, imnames1 = load_imgs("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/bonobo_volumes_large", (256,256,256), true)
images2, imnames2 = load_imgs("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/chimp_volumes_large", (256,256,256), true)
images = cat(images1, images2, dims=4)
images = Float32.(images)
images1 = nothing
images2 = nothing

landmark_data = read_landmarks("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/bonobo_volumes_large", 22, "@1")
landmark_data2 = read_landmarks("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/chimp_volumes_large", 22, "@1")
landmarks = cat(landmark_data, landmark_data2, dims=2)./10

X_train, X_test, y_train, y_test = regular_train_test_split_2d(images, landmarks)

X_train, y_train = (mirror_vol(X_train, y_train))

X_train = Float32.(X_train)

flip1, lm_flip1 = flip_3D(X_train, y_train)
flip2, lm_flip2 = flip_volume_front(flip1, lm_flip1)
flip3, lm_flip3 = flip_volume_side(flip2, lm_flip2)

rot, lm_rot = rotate_volumes(flip2, lm_flip2, 10)
rot2, lm_rot2 = rotate_volumes(X_train, y_train, 20)
rot3, lm_rot3 = rotate_volumes(flip3, lm_flip3, -15)

jit, lm_jit = jitter_3D(flip1, lm_flip1, 10)
jit2, lm_jit2 = jitter_3D(X_train, y_train, 10)
jit3, lm_jit3 = jitter_3D(flip3, lm_flip3, 10)
jit4, lm_jit4 = jitter_3D(rot2, lm_rot2, 10)

using ImageView
ImageView.imshow(X_test)

complete1 = depth_map_all_sides(X_train)
X_train=nothing
complete2 = depth_map_all_sides(flip1)
flip1=nothing
complete3 = depth_map_all_sides(flip2)
flip2=nothing
complete4 = depth_map_all_sides(flip3)
flip3=nothing
complete5 = depth_map_all_sides(rot)
rot=nothing
complete6 = depth_map_all_sides(rot2)
rot2=nothing
complete7 = depth_map_all_sides(rot3)
rot3=nothing
complete8 = depth_map_all_sides(jit)
jit=nothing
complete9 = depth_map_all_sides(jit2)
jit2=nothing
complete10 = depth_map_all_sides(jit3)
jit3=nothing
complete11 = depth_map_all_sides(jit4)

X_train = cat(complete1, complete2, complete3, complete4, complete5, complete6,
    complete7, complete8, complete9, complete10, complete11, dims=4)
y_train = cat(y_train, lm_flip1, lm_flip2, lm_flip3, lm_rot, lm_rot2, lm_rot3, lm_jit, lm_jit2, lm_jit3, lm_jit4, dims=2)

X_train = image_gradients(X_train)
# make depthmaps of the testing data
X_test = depth_map_all_sides(X_test)
X_test = image_gradients(X_test)

# define the cost function
cost(x, y) = sum((model(x)-y).^2)|>gpu

#define the model
model = Flux.mapleaves(cu, AutoLandmarking.vgg19)

# define the trainingrate and optimiser
opt = Flux.ADAM(0.000015)

# redefinition of the dropout function to work with testmode!() since this
# functionality is not working at the moment. This step only needs to be done
# if you wish to perform uncertainty estimation using the dropout method (AutoLM.response_distribution())
# after training.
using Random
function Flux.dropout(x, p; dims = :)
    q = 1 - p
    y = rand!(similar(x, Flux._dropout_shape(x, dims)))
    y .= Flux._dropout_kernel.(y, p, q)
    x .* y
end

import Zygote
Zygote.@adjoint function Flux.dropout(x, p; dims = :)
   q = 1 - p
   y = rand!(similar(x, Flux._dropout_shape(x, dims)))
   y .= Flux._dropout_kernel.(y, p, q)
   return x .* y, Δ -> (Δ .* y, nothing)
end
Zygote.refresh()

function run_model(modell, X, y)
  train_data = Flux.mapleaves(cu, X)
  train_labels = Flux.mapleaves(cu, y)
  dataset = Flux.Data.DataLoader(train_data, train_labels, batchsize = 4, shuffle=true)
  Flux.train!(cost, params(model), dataset, opt)
  testmode!(model)
  cosima = cost_whole_data_2D(train_data, train_labels, cost)
  testmode!(model, false)
  return cosima
end

# create lists to store metrics
accs = []
costs = []

# run the training for 300 epochs. batches of 128 will be loaded onto the GPU,
# which will be further subdevided into minibatches of 4 (as defined in AutoML.run_model).
for i in 1:150
  costr = 0
  for j in 1:128:size(X_train, 4)
    if j+127 > size(X_train, 4)
      train_data = X_train[:,:,:,j:end]
      train_labels = y_train[:,j:end]
    else
      train_data = X_train[:,:,:,j:j+127]
      train_labels = y_train[:,j:j+127]
    end
    cosima = run_model(model, train_data, train_labels)
    costr = costr + cosima
    train_data = nothing
    train_labels = nothing
    # GC.gc()
  end
  push!(costs, costr)
  println("epoch ", i, " finished")
  println(costr)
  testmode!(model)
  acc2, max1, min1 = AutoLandmarking.avg_accuracy_per_point(model, gpu(X_test[:,:,:,:]), gpu(y_test[:,:]), 3)
  testmode!(model, false)
  println("median deviation per point on testing dataset: ", acc2, "with maximum", max1)
  println("median sum is ", sum(acc2))
  push!(accs, sum(acc2))
  if sum(acc2)[1]<23.0
    break
  end
end

#solve issue with arrays inside vector
aks = []
for i in accs
  push!(aks, i[1])
end

# check what the best accuracy is
minimum(aks)

#plot the development of cost and accuracy
import Plots
Plots.plot(aks, legend=:topright, label="sum of deviations per point", color= :red, xlabel="epochs", ylabel = "acuracy")
plt = Plots.plot!(Plots.twinx(), costs, label="training loss", legend=:topleft, ylabel = "loss")
Plots.png(plt, "C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/mandibles_training.png")

# Data analysis mandibles ==================
response = cpu(predict_set(gpu(X_test), model))
X_train2, X_test2, y_train2, y_test2 = regular_train_test_split_2d(images, landmarks)
on_vol = AutoLandmarking.landmark_to_surface(X_test2, response, 6)
on_vol_no4 = on_vol[:, cat(1:3,5:19,dims=1)]
response_no4 = response[:, cat(1:3,5:19,dims=1)]

function print_costs(y, resp)
  costs=[]
  for i in 1:size(resp,2)
    dev = sum((y[:,i] .- resp[:,i]) .^2)
    push!(costs, dev)
    println(i, "    ", dev)
  end
  return costs
end

costs_resp = print_costs(y_test_no4, response_no4)

# puting on volume makes it worse

three_d_train = AutoLandmarking.to_3d_array(y_train2)
three_d_test = AutoLandmarking.to_3d_array(response_no4)

aligned_train = AutoLandmarking.align_all(three_d_train)
aligned_test = AutoLandmarking.align_all(three_d_test)

mean_train, stdevs = mean_shape(aligned_train)

distances = proc_distance(mean_train, aligned_test, stdevs)
dists = []
for i in 1:18
  push!(dists, distances[i])
end

for i in 1:18
  println("ind $i cost is ", costs_resp[i], "  and dist is ", dists[i])
end

uncertainties = []
for i in 1:18
  means, stddev = AutoLandmarking.response_distribution(model, gpu(X_test_no4[:,:,:,i:i]), gpu(y_test_no4[:,i:i]), 100)
  println("$i has uncevertaty", sum(stddev), "  and cost  ", costs_resp[i])
  push!(uncertainties, sum(stddev))
end
# 6 are over 40

dists_conf = dists[findall(x->x<4,dists)]
resp_conf = response_no4[:,findall(x->x<4,dists)]
volumes_conf = X_test2[:,:,:,findall(x->x<4,dists)]
y_test_conf = y_test_no4[:,findall(x->x<4,dists)]
X_test_conf = X_test_no4[:,:,:,findall(x->x<4,dists)]
AutoLandmarking.change_values!(y_test_conf, 26.0, 0, >)
ImageView.imshow(X_test_conf)

is_bonobo = zeros(97,1)
is_bonobo[1:37,1] .= 1
imnames = cat(imnames1, imnames2, dims=1)
test_names = []
train_names = []
is_bonobo_train = []
is_bonobo_test = []
for i in 1:97
  if i % 5 == 0
    push!(test_names, imnames[i])
    push!(is_bonobo_test, is_bonobo[i])
  else
    push!(train_names, imnames[i])
    push!(is_bonobo_train, is_bonobo[i])
  end
end

is_bonobo_test_no4 = is_bonobo_test[cat(1:3,5:19,dims=1)]
is_bonobo_conf = is_bonobo_test_no4[findall(x->x<41,dists)]
test_names = test_names[cat(1:3,5:19,dims=1)]
names_conf = test_names[findall(x->x<41,dists)]
AutoLandmarking.array_to_lm_file("C:/Users/immanueldiener/Desktop/Master/worst_pred.landmarkAscii", response_no4[:,18])

import ImageView
using DelimitedFiles
using Statistics

y_test_no4 = y_test[:,cat(1:3,5:19,dims=1)]
X_test_no4 = X_test[:,:,:,cat(1:3,5:19,dims=1)]
ImageView.imshow(X_test_no4)

ImageView.imshow(depth_map_all_sides(X_test_no4))
test_names[14]
all_devs = zeros(14, 4)
devs_per_point = zeros(14*22, 2)

AutoLandmarking.change_values!(y_test_no4, 26.0, 0.0, >)

testmode!(model)
for i in 1:14
  acc2, max1, min1 = AutoLandmarking.avg_accuracy_per_point(model, gpu(X_test_conf[:,:,:,i:i]), gpu(y_test_conf[:,i:i]), 3)
  all_devs[i,1] = sum(acc2)
  all_devs[i,2] = dists_conf[i]
  all_devs[i,3] = i
  if i in findall(x->x<41,dists)
    all_devs[i,4] = 1
  end
end
minimum(all_devs[:,1])
all_devs
using Distances
for i in 1:14
  for j in 1:22
    devs_per_point[(i-1)*22+j, 1] = euclidean(resp_conf[j*3-2:j*3, i], y_test_conf[j*3-2:j*3, i])
    devs_per_point[(i-1)*22+j, 2] = j
    if j == 2
      println(i, "   ", euclidean(resp_conf[j*3-2:j*3, i], y_test_conf[j*3-2:j*3, i]))
    end
  end
end

devs_per_point

AutoLandmarking.change_values!(y_test, 25.8, 25.6, >)
test_names[18]

out = hcat(all_devs, names_conf)

writedlm("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/devs_per_ind.txt", out)
writedlm("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/devs_per_point.txt", devs_per_point)

# timing GPU : 12000/min, CPU: 110/min
@time begin
  predict_set(X_train[:,:,:,1:50], model_cpu)
end

# export of coordinates
out_test = hcat(names_conf, is_bonobo_conf, resp_conf')
out_train = hcat(train_names, is_bonobo_train, y_train')
out_all = vcat(out_train, out_test)
is_train = zeros(92, 1)
is_train[1:78] .= 1
output_all = hcat(out_all[:,1:2], is_train, out_all[:,3:end])

writedlm("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/coords_pred.txt", out_test)
writedlm("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/coords_train.txt", out_train)
writedlm("C:/Users/immanueldiener/Desktop/Master/master_data/mandibles/coords_all.txt", output_all)
# ===========================================

using BSON
testmode!(model)
model_cpu = cpu(model)
BSON.@save "C:/Users/immanueldiener/Desktop/Master/master_data/model_mandibles.bson" model_cpu
BSON.@load "C:/Users/immanueldiener/Desktop/Master/master_data/model_mandibles.bson" model_cpu
model = gpu(model_cpu)

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
"""
    galerkin_matrix(
        B::AbstractBSplineBasis,
        [deriv = (Derivative(0), Derivative(0))],
        [MatrixType = BandedMatrix{Float64}],
    )

Compute Galerkin mass or stiffness matrix, as well as more general variants
of these.

# Extended help

The Galerkin mass matrix is defined as

```math
M_{ij} = ⟨ ϕ_i, ϕ_j ⟩ \\quad \\text{for} \\quad
i ∈ [1, N] \\text{ and } j ∈ [1, N],
```

where ``ϕ_i(x)`` is the ``i``-th basis function and `N = length(B)` is the
number of functions in the basis `B`.
Here, ``⟨f, g⟩`` is the [``L^2`` inner
product](https://en.wikipedia.org/wiki/Square-integrable_function#Properties)
between functions ``f`` and ``g``.

Since products of B-splines are themselves piecewise polynomials, integrals can
be computed exactly using [Gaussian quadrature
rules](https://en.wikipedia.org/wiki/Gaussian_quadrature).
To do this, we use Gauss--Legendre quadratures via the
[FastGaussQuadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
package.

## Matrix layout and types

The mass matrix is banded with ``2k - 1`` bands.
Moreover, the matrix is symmetric and positive definite, and only ``k`` bands are
needed to fully describe the matrix.
Hence, a
[`Hermitian`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.Hermitian)
view of an underlying matrix is returned.

By default, the underlying matrix holding the data is a `BandedMatrix` that
defines the upper part of the symmetric matrix.
Other types of container are also supported, including regular sparse matrices
(`SparseMatrixCSC`) and dense arrays (`Matrix`).
See [`collocation_matrix`](@ref) for a discussion on matrix types.

## Derivatives of basis functions

Galerkin matrices associated to the derivatives of basis functions may be
constructed using the optional `deriv` parameter.
For instance, if `deriv = (Derivative(0), Derivative(2))`, the matrix
``⟨ ϕ_i, ϕ_j'' ⟩`` is constructed, where primes denote derivatives.
Note that, if the derivative orders are different, the resulting matrix is not
symmetric, and a `Hermitian` view is not returned in those cases.

## Combining different bases

More generally, it is possible to compute matrices of the form
``⟨ ψ_i^{(n)}, ϕ_j^{(m)} ⟩``, where `n` and `m` are derivative orders, and ``ψ_i``
and ``ϕ_j`` belong to two *different* (but related) bases `B₁` and `B₂`.
For this, instead of the `B` parameter, one must pass a tuple of bases
`(B₁, B₂)`.
The restriction is that the bases must have the same parent B-spline basis.
That is, they must share the same set of B-spline knots and be of equal
polynomial order.

Note that, if both bases are different, the matrix will not be symmetric, and
will not even be square if the bases have different lengths.

In practice, this feature may be used to combine a B-spline basis `B`, with a
recombined basis `R` generated from `B` (see [Basis recombination](@ref basis-recombination-api)).
"""
function galerkin_matrix(
        Bs::NTuple{2,AbstractBSplineBasis},
        deriv::DerivativeCombination{2} = Derivative.((0, 0)),
        ::Type{M} = BandedMatrix{Float64},
    ) where {M <: AbstractMatrix}
    B1, B2 = Bs
    symmetry = M <: Hermitian
    if symmetry && (deriv[1] !== deriv[2] || B1 !== B2)
        throw(ArgumentError(
            "input matrix type incorrectly assumes symmetry"))
    end
    A = allocate_galerkin_matrix(M, Bs)
    galerkin_matrix!(A, Bs, deriv)
end

function galerkin_matrix(
        B::AbstractBSplineBasis,
        deriv::DerivativeCombination{2} = Derivative.((0, 0)),
        ::Type{Min} = BandedMatrix{Float64},
    ) where {Min <: AbstractMatrix}
    symmetry = deriv[1] === deriv[2]
    T = eltype(Min)
    M = symmetry ? Hermitian{T,Min} : Min
    galerkin_matrix((B, B), deriv, M)
end

galerkin_matrix(B, ::Type{M}) where {M <: AbstractMatrix} =
    galerkin_matrix(B, Derivative.((0, 0)), M)

function _check_bases(Bs::Tuple{Vararg{AbstractBSplineBasis}})
    ps = parent.(Bs)
    if any(ps .!== ps[1])
        throw(ArgumentError(
            "all bases must share the same parent B-spline basis"
        ))
    end
    nothing
end

allocate_galerkin_matrix(::Type{Hermitian{T,M}},
                         Bs) where {T, M <: AbstractMatrix} =
    Hermitian(allocate_galerkin_matrix(M, Bs, true))

allocate_galerkin_matrix(::Type{M}, Bs) where {M} =
    allocate_galerkin_matrix(M, Bs, false)

allocate_galerkin_matrix(::Type{M}, Bs, _symmetry) where {M <: AbstractMatrix} =
    M(undef, length.(Bs)...)

allocate_galerkin_matrix(::Type{SparseMatrixCSC{T}}, Bs, _symmetry) where {T} =
    spzeros(T, length.(Bs)...)

function allocate_galerkin_matrix(::Type{M}, Bs, symmetry) where {M <: BandedMatrix}
    _check_bases(Bs)
    B1, B2 = Bs
    k = order(B1)
    @assert k == order(B2)  # verified in _check_bases
    # When the bases are the same, the upper/lower bandwidths are Nb = k - 1
    # (total = 2k - 1 bands).
    # If the matrix is symmetric, we only store the upper band.
    Nb = k - 1
    N1, N2 = length.(Bs)
    bands = if symmetry
        (0, Nb)
    else
        # If the bases have different number of constraints (BCs) on the left
        # boundary (⇔ δl ≠ 0), then the bands are shifted.
        δl, δr = num_constraints(B2) .- num_constraints(B1)
        @assert N1 == N2 + δl + δr
        (Nb + δl, Nb - δl)
    end
    M(undef, (N1, N2), bands)
end

"""
    galerkin_matrix!(A::AbstractMatrix, B::AbstractBSplineBasis,
                     deriv = (Derivative(0), Derivative(0)))

Fill preallocated Galerkin matrix.

The matrix may be a `Hermitian` view, in which case only one half of the matrix
will be filled. Note that, for the matrix to be symmetric, both derivative orders
in `deriv` must be the same.

More generally, it is possible to combine different functional bases by passing
a tuple of `AbstractBSplineBasis` as `B`.

See [`galerkin_matrix`](@ref) for details.
"""
function galerkin_matrix! end

galerkin_matrix!(M, B::AbstractBSplineBasis, args...) =
    galerkin_matrix!(M, (B, B), args...)

function galerkin_matrix!(
        S::AbstractMatrix, Bs::Tuple{Vararg{<:AbstractBSplineBasis,2}},
        deriv = Derivative.((0, 0)),
    )
    _check_bases(Bs)
    B1, B2 = Bs
    same_ij = B1 == B2 && deriv[1] == deriv[2]
    T = eltype(S)

    if size(S) != length.(Bs)
        throw(DimensionMismatch("wrong dimensions of Galerkin matrix"))
    end

    # Orders and knots are assumed to be the same (see _check_bases).
    k = order(B1)
    ts = knots(B1)
    @assert k == order(B2) && ts === knots(B2)

    # Quadrature information (nodes, weights).
    quadx, quadw = _quadrature_prod(Val(2k - 2))
    @assert length(quadx) == k  # we need k quadrature points per knot segment

    A = if S isa Hermitian
        deriv[1] === deriv[2] ||
            throw(ArgumentError("matrix will not be symmetric with deriv = $deriv"))
        B1 === B2 ||
            throw(ArgumentError("matrix will not be symmetric if bases are different"))
        fill_upper = S.uplo === 'U'
        fill_lower = S.uplo === 'L'
        parent(S)
    else
        fill_upper = true
        fill_lower = true
        S
    end

    fill!(S, 0)
    nlast = last(eachindex(ts))
    ioff = first(num_constraints(B1))
    joff = first(num_constraints(B2))

    # We loop over all knot segments Ω[n] = (ts[n], ts[n + 1]).
    # We integrate all supported B-spline products B1[i] * B2[j] over this
    # segment, adding the result to A[i, j].
    @inbounds for n in eachindex(ts)
        n == nlast && break
        tn, tn1 = ts[n], ts[n + 1]
        tn1 == tn && continue  # interval of length = 0

        metric = QuadratureMetric(tn, tn1)

        # Unnormalise quadrature nodes, such that xs ∈ [tn, tn1]
        xs = metric .* quadx
        # @assert all(x -> tn ≤ x ≤ tn1, xs)

        for (x, w) ∈ zip(xs, quadw)
            ilast = n - ioff
            jlast = n - joff
            _, bis = evaluate_all(B1, x, deriv[1], T; ileft = ilast)
            _, bjs = same_ij ?
                (ilast, bis) : evaluate_all(B2, x, deriv[2], T; ileft = jlast)
            y = metric.α * w
            for (δj, bj) ∈ pairs(bjs), (δi, bi) ∈ pairs(bis)
                i = ilast + 1 - δi
                j = jlast + 1 - δj
                if !fill_upper && i < j
                    continue
                elseif !fill_lower && i > j
                    continue
                elseif i ∉ axes(A, 1)  # this can be true for recombined bases
                    @assert iszero(bi)
                    continue
                elseif j ∉ axes(A, 2)
                    @assert iszero(bj)
                    continue
                end
                @inbounds A[i, j] += y * bi * bj
            end
        end
    end

    S
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
"""
    BoseHubbardReal1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0)

Implements a two-component one-dimensional Bose Hubbard chain in real space.

```math
\\hat{H} = \\hat{H}_a + \\hat{H}_b + V\\sum_{i} n_{a_i}n_{b_i}
```

# Arguments

* `address`: the starting address, defines number of particles and sites.
* `ua`: the on-site interaction parameter parameter for Hamiltonian a.
* `ub`: the on-site interaction parameter parameter for Hamiltonian b.
* `ta`: the hopping strength for Hamiltonian a.
* `tb`: the hopping strength for Hamiltonian b.
* `v`: the inter-species interaction parameter V.

# See also

* [`HubbardReal1D`](@ref)
* [`BoseHubbardMom1D2C`](@ref)

"""
struct BoseHubbardReal1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
    ha::HA
    hb::HB
end

function BoseHubbardReal1D2C(add::BoseFS2C; ua=1.0,ub=1.0,ta=1.0,tb=1.0,v=1.0)
    ha = HubbardReal1D(add.bsa; u=ua, t=ta)
    hb = HubbardReal1D(add.bsb; u=ub, t=tb)
    T = promote_type(eltype(ha), eltype(hb))
    return BoseHubbardReal1D2C{T,typeof(ha),typeof(hb),v}(ha, hb)
end

function Base.show(io::IO, h::BoseHubbardReal1D2C)
    addr = starting_address(h)
    ua = h.ha.u
    ub = h.hb.u
    ta = h.ha.t
    tb = h.hb.t
    v = h.v
    print(io, "BoseHubbardReal1D2C($addr; ua=$ua, ub=$ub, ta=$ta, tb=$tb, v=$v)")
end

function starting_address(h::BoseHubbardReal1D2C)
    return BoseFS2C(starting_address(h.ha), starting_address(h.hb))
end

LOStructure(::Type{<:BoseHubbardReal1D2C{<:Real}}) = IsHermitian()

Base.getproperty(h::BoseHubbardReal1D2C, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::BoseHubbardReal1D2C, ::Val{:ha}) = getfield(h, :ha)
Base.getproperty(h::BoseHubbardReal1D2C, ::Val{:hb}) = getfield(h, :hb)
Base.getproperty(h::BoseHubbardReal1D2C{<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V

# number of excitations that can be made
function num_offdiagonals(ham::BoseHubbardReal1D2C, add)
    return 2*(num_occupied_modes(add.bsa) + num_occupied_modes(add.bsb))
end

"""
    bose_hubbard_2c_interaction(::BoseFS2C)

Compute the interaction between the two components.
"""
function bose_hubbard_2c_interaction(add::BoseFS2C)
    c1 = onr(add.bsa)
    c2 = onr(add.bsb)
    interaction = zero(eltype(c1))
    for site = 1:length(c1)
        if !iszero(c2[site])
            interaction += c2[site] * c1[site]
        end
    end
    return interaction
end

function diagonal_element(ham::BoseHubbardReal1D2C, address::BoseFS2C)
    return (
        diagonal_element(ham.ha, address.bsa) +
        diagonal_element(ham.hb, address.bsb) +
        ham.v * bose_hubbard_2c_interaction(address)
    )
end

function get_offdiagonal(ham::BoseHubbardReal1D2C, add, chosen)
    nhops = num_offdiagonals(ham,add)
    nhops_a = 2 * num_occupied_modes(add.bsa)
    if chosen ≤ nhops_a
        naddress_from_bsa, onproduct = hopnextneighbour(add.bsa, chosen)
        elem = - ham.ha.t*sqrt(onproduct)
        return BoseFS2C(naddress_from_bsa,add.bsb), elem
    else
        chosen -= nhops_a
        naddress_from_bsb, onproduct = hopnextneighbour(add.bsb, chosen)
        elem = -ham.hb.t * sqrt(onproduct)
        return BoseFS2C(add.bsa,naddress_from_bsb), elem
    end
    # return new address and matrix element
end

struct OffdiagonalsBoseReal1D2C{
    A<:BoseFS2C,T,H<:TwoComponentHamiltonian{T}
} <: AbstractOffdiagonals{A,T}
    hamiltonian::H
    address::A
    length::Int
    num_hops_a::Int
end

function offdiagonals(h::BoseHubbardReal1D2C, a::BoseFS2C)
    hops_a = num_offdiagonals(h.ha, a.bsa)
    hops_b = num_offdiagonals(h.hb, a.bsb)
    length = hops_a + hops_b

    return OffdiagonalsBoseReal1D2C(h, a, length, hops_a)
end

function Base.getindex(s::OffdiagonalsBoseReal1D2C{A,T}, i)::Tuple{A,T} where {A,T}
    @boundscheck 1 ≤ i ≤ s.length || throw(BoundsError(s, i))
    if i ≤ s.num_hops_a
        new_a, matrix_element = get_offdiagonal(s.hamiltonian.ha, s.address.bsa, i)
        new_add = A(new_a, s.address.bsb)
    else
        i -= s.num_hops_a
        new_b, matrix_element = get_offdiagonal(s.hamiltonian.hb, s.address.bsb, i)
        new_add = A(s.address.bsa, new_b)
    end
    return new_add, matrix_element
end

Base.size(s::OffdiagonalsBoseReal1D2C) = (s.length,)


text:
# ---
# title: 1493. Longest Subarray of 1's After Deleting One Element
# id: problem1493
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/longest-subarray-of-1s-after-deleting-one-element/description/>
# hidden: true
# ---
#
# Given a binary array `nums`, you should delete one element from it.
#
# Return the size of the longest non-empty subarray containing only 1's in the
# resulting array.
#
# Return 0 if there is no such subarray.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,0,1]
#     Output: 3
#     Explanation: After deleting the number in position 2, [1,1,1] contains 3 numbers with value of 1's.
#
# **Example 2:**
#
#
#
#     Input: nums = [0,1,1,1,0,1,1,0,1]
#     Output: 5
#     Explanation: After deleting the number in position 4, [0,1,1,1,1,1,0,1] longest subarray with value of 1's is [1,1,1,1,1].
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1,1]
#     Output: 2
#     Explanation: You must delete one element.
#
# **Example 4:**
#
#
#
#     Input: nums = [1,1,0,0,1,1,1,0,1]
#     Output: 4
#
#
# **Example 5:**
#
#
#
#     Input: nums = [0,0,0]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `nums[i]` is either `0` or `1`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Statistics, Random

# version which generates shock internally
function SVmodel(θ, n, burnin)
    ϕ = θ[1]
    ρ = θ[2]
    σ = θ[3]
    hlag = 0.0
    ys = zeros(n)
    for t = 1:burnin+n
        # bound the log variance at a high value,
        # too bullet-proof things which use this
        h = min(ρ*hlag + σ*randn(),20.0)
        y = ϕ*exp(h/2.0)*randn()
        if t > burnin
            ys[t-burnin] = y
        end
        hlag = h
    end
    ys, ϕ*exp(hlag/2.0) # return the sample of returns, plus the final period volatility
end
# the dgp: simple discrete time stochastic volatility (SV) model
function SVmodel(θ, n, η, ϵ, savedata=false)
    ϕ = θ[1]
    ρ = θ[2]
    σ = θ[3]
    burnin = size(η,1) - n
    hlag = 0.0
    ys = zeros(n,1)
    for t = 1:burnin+n
        h = ρ*hlag + σ*η[t]
        σt = ϕ*exp(h/2.0)
        y = σt*ϵ[t]
        if t > burnin
            ys[t-burnin] = y
        end
        hlag = h
    end
    if savedata == true
        writedlm("svdata.txt", ys)
    end
    ys, ϕ*exp(hlag/2.0) # return the sample of returns, plus the final period volatility
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
function radixgroup!(x::Vector{T}, uint::UInt = UInt(1)) where T<:AbstractString
    @time hashindex = hash.(x)
    @time index = collect(1:length(x))
    @time grouptwo!(hashindex, index)
    @time x_permuted = @view x[index]
    @time if isgrouped(x_permuted, hashindex)
        return index
    else
        return radixgroup!(x, uint + 1)
    end
end

function isgrouped(x, hashx)
    x1 = x[1]
    hashx1 = hashx[1]
    for i = 2:length(x)
        @inbounds if x1 != x[i]
            if hashx1 == hashx[i]
                return false
            else
                x1 = x[i]
                hashx1 = hashx[i]
            end
        end
    end
    return true
end

# @time permi = radixgroup!(idstr)

function fastby(fn::Function, byvec::Vector{T}, valvec::Vector) where T <: AbstractString
    permi = radixgroup!(byvec)
    fastby_contiguous(fn, byvec[permi], valvec[permi])
end

function fastby_contiguous(fn::Function, byvec::Vector{T}, valvec::Vector{S}) where {T <: AbstractString, S}
    res = Dict{T, S}()
    l = length(byvec)

    j = 1
    lastby = byvec[1]
    for i = 2:l
        @inbounds byval = byvec[i]
        if byval != lastby
            viewvalvec = @view valvec[j:i-1]
            @inbounds res[lastby] = fn(viewvalvec)
            j = i
            @inbounds lastby = byvec[i]
        end
    end

    viewvalvec = @view valvec[j:l]
    @inbounds res[byvec[l]] = fn(viewvalvec)
    return res
end

@time fastby(sum, idstr, valvec)

using StatsBase
@time countmap(idstr)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

# XXX: trailsign keyword not implemented, not sure it's worth the effort.
# Possible strategy: hold the sign of number in a variable (you can use
# copysign(1, number)), set result[1] equal to dd and multiply the first
# non-zero element (use findfirst for that) by the sign of number.
function _sixty{T<:AbstractFloat}(number::T)
    dd = trunc(abs(number))
    mm = abs(60.0*number)
    ss = abs(3600.0*number)
    result = zeros(typeof(number), 3)
    result[1] = trunc(number)
    result[2] = trunc(mm - 60.0*dd)
    result[3] = ss - 3600.0*dd - 60.0*result[2]
    return result
end

"""
    sixty(number) -> [deg, min, sec]

### Purpose ###

Converts a decimal number to sexagesimal.

### Explanation ###

The reverse of `ten` function.

### Argument ###

* `number`: decimal number to be converted to sexagesimal.

### Output ###

An array of three `AbstractFloat`, that are the sexagesimal counterpart
(degrees, minutes, seconds) of `number`.

### Example ###

``` julia
sixty(-0.615)
# => 3-element Array{Float64,1}:
#     -0.0
#     36.0
#     54.0
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
sixty(n::Real) = _sixty(float(n))


text:
using HorizonSideRobots

function stairs_marking_with_obstacles(r::Robot)
    back_track = go_to_bottom_left_corner_return_path(r)

    putmarker!(r)
    while isborder(r, Ost) == false
        move!(r, Ost)
        putmarker!(r)
    end
    width = show_path_after_moving(r, West)

    while isborder(r, Nord) == false
        width_temp = width
        move!(r, Nord)
        while width_temp > 0
            putmarker!(r)
            width_temp -= move_through_obstacle_and_return_step_count(r, Ost)
        end

        step_count = 1
        while step_count !=0
            step_count = move_through_obstacle_and_return_step_count(r, West)
        end

        width -= 1
    end

    go_to_bottom_left_corner_return_path(r)
    back_by_path_o_sides(r, back_track)
end

function go_to_bottom_left_corner_return_path(r::Robot)::Array
    path = []
    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        if isborder(r, Sud) == false
            move!(r, Sud)
            push!(path, Nord)
        end
        if isborder(r, West) == false
            move!(r, West)
            push!(path, Ost)
        end
    end
    return path
end

function back_by_path_o_sides(r, path::Array)
    n = length(path)
    while n > 0
        move!(r, path[n])
        n = n - 1
    end
end

function move_through_obstacle_and_return_step_count(r::Robot, side::HorizonSide)
    n = 0
    side_side = HorizonSide(mod(Int(side)+1,4)) # Следующее направление против часовой стрелки
    while isborder(r, side) && (isborder(r, side_side) == false) # Следующее против часовой стрелки
        move!(r, side_side)
        n += 1
    end
    step_count = 0
    if isborder(r, side) == false
        move!(r, side)
        step_count += 1
    end
    if n != 0
        while isborder(r, HorizonSide(mod(Int(side_side) + 2, 4)))
            move!(r,side)
            step_count +=1
        end
        for i in 1:n
            move!(r, HorizonSide(mod(Int(side_side) + 2, 4)))
        end
    end
    return step_count
end

function show_path_after_moving(r::Robot, side::HorizonSide)
    n = 0
    while isborder(r, side)==false
        move!(r, side)
        n = n + 1
    end
    return n
end


text:
using Revise
using Flipping
using Guilia
using MixedModels
using Distributions
##
Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Test"
Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Dev_raw_data"
pokes, streaks, DataIndex = create_exp_dataframes(Directory)
###
Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Dev_raw_data/Results/"
file = joinpath(Directory,"streaksResults.csv")
streaks = Guilia.carica(file)
exclude_mice = ["CN21", "CN22","CN23","CN24","CN41","CN42","CN43"]# CNO GROUP WITHOUT CNO THE FIRST DAy
streaks = @filter streaks !in(:MouseID, exclude_mice)
Ncs = ["NC"* string(i) for i in collect(1:14)]
Ccs = ["CC"* string(i) for i in [3,4,5,9,10,14,15,16]]
Rcs = ["RC"* string(i) for i in [1,2,3,4]]
Nbs = ["NB"* string(i) for i in [21,22,23,24,25,41,42,43,44,45]]
Bns = ["BN"* string(i) for i in [21,22,23,24,41,42,43,44]]
Its = ["IT"* string(i) for i in [21,22,23,24,41,42,43,44]]
age_exp = vcat(Nbs,Bns,Its)
youngs = ["BN21", "BN22", "BN41", "BN42", "IT21", "IT22","IT41", "IT42", "NB21", "NB22","NB23","NB41","NB42"]
dev_below15_day1 = ["BN21","BN23","BN41","BN42","BN44","IT23","NB25","NB43","NB45"]
proj_below15_day1 = ["CC9","NC10","NC12","NC14","RC3","RC8"]
below_thrs = vcat(dev_below15_day1,proj_below15_day1)
cnos_animals = vcat(Ncs,Ccs,Rcs)
wt =  append!(["NC"* string(i) for i in [1,2,3,4,9,10,11,12]],age_exp)
"BN21" in youngs
streaks = @apply streaks begin
    @transform {Age = :MouseID in youngs ? "Young" : "Adult"}
    @transform {Exp_type = :MouseID in age_exp ? "Development" : "Projections"}
    @transform {Treatment = :MouseID in cnos ? "CNO" : "SAL"}
    @transform {Gen = :MouseID in wt ? "WT" : "HET"}
    @transform {Low_performance = :MouseID in below_thrs}
    @transform {Combo = :Gen*"_"*:Treatment}
end

union(column(streaks,:MouseID))
union(column(streaks,:Age))
filetosave = joinpath(Directory,"Full_info_streaks.csv")
CSV.write(filetosave,streaks)

cno_only = @apply streaks begin
    @filter !in(:MouseID,age_exp)
    @filter :Exp_Day == 1
end

check = @filter cno_only :Treatment == "CNO"
filetosave = joinpath(Directory,"All_dreadds_streaks.csv")
CSV.write(filetosave,streaks)
#######
data = launch()
t = data[][:data][]
df = DataFrame(t)

categorical!(df,:Treatment)
categorical!(df,:Gen)
##
function Likelyhood_Ratio_test(simple,full)
    degrees = dof(full) - dof(simple)
    ccdf(Distributions.Chisq(degrees), deviance(simple) - deviance(full))
end

function AIC_test(candidate, simpler)
    exp((aic(candidate) - aic(simpler))/2)
end
function AICc(model)
    aic(model) + ((2*(dof(model)^2) + 2*dof(model))/(nobs(model) - dof(model) - 1))
end
function AICc_test(candidate, simpler)
    exp((AICc(candidate) - AICc(simpler))/2)
end
##
No_interaction = fit(LinearMixedModel, @formula( AfterLast ~ 1 + Treatment + Gen + (1|MouseID)), df)
Interaction = fit(LinearMixedModel, @formula( AfterLast ~ 1 + Treatment + Gen + Treatment * Gen + (1|MouseID)), df)

AIC_test(protocol_effect,zero_effects)
AICc_test(protocol_effect,zero_effects)
Likelyhood_Ratio_test(zero_effects,protocol_effect)


text:
using DataFrames
using Plots
using Query
using StatPlots
using MarketTechnicals
using GLMNet
using GLM

cases = readtable("SCDB_2017_01_caseCentered_Citation.csv")
justices = readtable("SCDB_2017_01_justiceCentered_Citation.csv")

function analyze()
     years = zeros(Int, size(cases)[1])
     for k = 1:length(years)
          years[k] = parse(Int, cases[:docketId][k][1:4])
     end
     cases[:year] = years

     plot(title = "Number of Laws Declared Unconstitutional", xlab = "Year", ylab = "Frequency",  bgcolor=RGBA(0.0,0.0,0.0,0))
     topics = ["No declaration of unconstitutionality",
               "Act of Congress",
               "State or Territorial Law",
               "Municipal or Local Ordinance"]
     colors = distinguishable_colors(10)[4:8]
     for topic = 2:4
          uncon = cases |> @filter(_.declarationUncon == topic) |> DataFrame
          freq = uncon |> @groupby(_.year) |> @map({year=_.key, num_cases=length(_)}) |> DataFrame
          plot!(freq[:year], freq[:num_cases], label = topics[topic], linecolor = colors[topic])

          freq[:feat1] = freq[:year] - mean(freq[:year])
          freq[:feat2] = freq[:feat1].^2
          freq[:feat3] = freq[:feat1].^3
          freq[:val1] = freq[:num_cases] + 0.0
          model = lm(@formula(val1 ~ feat1 + feat2 + feat3), freq)
          plot!(freq[:year], GLM.predict(model), label = "", linestyle = :dash, linecolor = colors[topic])
     end
     png("Unconstitutional")
end

function topics()
     uncon = cases |> @filter(_.declarationUncon == 3) |> DataFrame
     labels = ["Criminal", "Civil", "1st", "Due", "Priv.", "Lawyer", "Union", "Econ.", "Jud.", "Fed.", "State", "Tax", "Misc", "Action"]

     freqs = countmap(uncon[:issueArea])
     sorted = []
     for k = 1:length(labels)
          append!(sorted, get(freqs, k, 0))
     end
     bar(labels, sorted, label = "", xlab = "Topics", ylab = "Number of Cases", title = "Topics in Unconstitutional Cases", bgcolor=RGBA(0.0,0.0,0.0,0))
     png("hist1")

     freqs = countmap(cases[:issueArea])
     sorted = []
     for k = 1:length(labels)
          append!(sorted, get(freqs, k, 0))
     end
     bar(labels, sorted, label = "", xlab = "Topics", ylab = "Number of Cases", title = "Topics in All Cases",  bgcolor=RGBA(0.0,0.0,0.0,0))
     png("hist2")
end

function rankings()
     df = DataFrame()
     df[:caseName] = cases[:caseName]
     df[:declarationUncon] = cases[:declarationUncon]

     caseNames = Dict{String, Integer}()
     for i = 1:size(cases)[1]
          name = cases[:caseName][i]
          caseNames[name] = i
     end

     freq = countmap(justices[:justiceName])
     for (k, v) in freq
          df[Symbol(k)] = -1
     end

     for i = 1:size(justices)[1]
          name = justices[:justiceName][i]
          case = justices[:caseName][i]
          vote = justices[:vote][i]
          if isna(vote)
               vote = 0
          end
          df[Symbol(name)][caseNames[case]] = vote
     end

     uncon = df |> @filter(_.declarationUncon == 3) |> DataFrame
     freq = countmap(justices[:justiceName])
     names = []
     for (k, v) in freq
          votes = countmap(uncon[Symbol(k)])
          yes = get(votes, 1, 0) + get(votes, 3, 0) + get(votes, 4, 0)
          no = get(votes, 2, 0) + get(votes, 6, 0) + get(votes, 7, 0)
          rat = yes/(yes+no)
          if rat < 0.5
               println(k, " ", yes, " ", no, " ", rat)
          end
     end
end

analyze()


text:
# native target for CPU execution

## target

export NativeCompilerTarget

Base.@kwdef struct NativeCompilerTarget <: AbstractCompilerTarget
    cpu::String=(LLVM.version() < v"8") ? "" : unsafe_string(LLVM.API.LLVMGetHostCPUName())
    features::String=(LLVM.version() < v"8") ? "" : unsafe_string(LLVM.API.LLVMGetHostCPUFeatures())
    always_inline::Bool=false # will mark the job function as always inline
end

llvm_triple(::NativeCompilerTarget) = Sys.MACHINE

function llvm_machine(target::NativeCompilerTarget)
    triple = llvm_triple(target)

    t = Target(triple=triple)

    tm = TargetMachine(t, triple, target.cpu, target.features)
    asm_verbosity!(tm, true)

    return tm
end

function process_entry!(job::CompilerJob{NativeCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
    if job.target.always_inline
        push!(function_attributes(entry), EnumAttribute("alwaysinline", 0, context(mod)))
    end
    invoke(process_entry!, Tuple{CompilerJob, LLVM.Module, LLVM.Function}, job, mod, entry)
end

## job

runtime_slug(job::CompilerJob{NativeCompilerTarget}) = "native_$(job.target.cpu)-$(hash(job.target.features))"


text:
# Static analysis of clamped twisted beam

# Source code: [`twisted_beam_tut.jl`](twisted_beam_tut.jl)

# ## Description

# The initially twisted cantilever beam is one of the standard test
# problems for verifying the finite-element accuracy [1]. The beam is
# clamped at one end and loaded either with unit in-plane or
# unit out-of-plane force at the other. The centroidal axis of the beam is
# straight at the undeformed  configuration, while its cross-sections are
# twisted about the centroidal axis from 0 at the clamped end to pi/2 at
# the free end.

# Reference deflection in the direction of the applied force.
# | Cross section thickness       | Loading in the Z direction | Loading in the Y direction |
# | -------   |     -------  |  ------- |
# | t = 0.32 | 0.005425   |    0.001753  |
# | t = 0.0032 |  0.005256  |  0.001294    |

# References:

# [1] MacNeal,  R. H., and R. L. Harder, “A Proposed Standard Set of Problems to
# Test Finite Element Accuracy,” Finite Elements in Analysis Design, vol. 11, pp.
# 3–20, 1985.

# [2] Simo,  J. C., D. D. Fox, and M. S. Rifai, “On a Stress Resultant Geometrically Exact Shell Model. Part II: The Linear Theory; Computational Aspects,” Computational Methods in Applied Mechanical Engineering, vol. 73, pp. 53–92, 1989.

# [3] Zupan D, Saje M (2004) On "A proposed standard set of problems to test
# finite element accuracy": the twisted beam. Finite Elements in Analysis
# and Design 40: 1445-1451.

# ## Goals

# # - Introduce definition of the shell model.
# # - Calculate the discrete model quantities and solve the static equilibrium problem.
# # - Demonstrate visualization of the resultant section forces and moments.

# ## Definition of the basic inputs

# The finite element code relies on the basic functionality implemented in this
# package.

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

# The inputs are defined in consistent units.
# The elastic properties are:

E = 0.29e8;
nu = 0.22;

# The material is elastic isotropic.
mater = MatDeforElastIso(DeforModelRed3D, E, nu)

# Here are the cross-section width and the cantilevered length of the beam.

W = 1.1;
L = 12.0;

# ## Reference Solutions

# The reference solutions are defined by these tuples of values: thickness,
# magnitude of the force, direction of the force, and the reference deflection
# along the force.

params_thicker_dir_3 = (t = 0.32, force = 1.0, dir = 3, uex = 0.005424534868469);
params_thicker_dir_2 = (t = 0.32, force = 1.0, dir = 2, uex = 0.001753248285256);

params_thinner_dir_3 = (t = 0.0032, force = 1.0e-6, dir = 3, uex = 0.005256);
params_thinner_dir_2 = (t = 0.0032, force = 1.0e-6, dir = 2, uex = 0.001294);

# Now select one particular simulation. Here we go with the thicker shell and
# the direction of the load 3:

params = params_thicker_dir_3

# ## Mesh generation

# The mesh is initially generated for a rectangular 2d domain. The element size
# can be controlled with these two variables:

nL = 48;
nW = 8;
nL, nW = 8, 4

fens, fes = T3block(L, W, nL, nW, :a);

# The 2d mesh is expanded into a three dimensional domain, and the
# locations of the nodes are tweaked to produce the pre twisted shape.
fens.xyz = xyz3(fens)
for i in 1:count(fens)
    a = fens.xyz[i, 1] / L * (pi / 2)
    y = fens.xyz[i, 2] - (W / 2)
    z = fens.xyz[i, 3]
    fens.xyz[i, :] = [fens.xyz[i, 1], y * cos(a) - z * sin(a), y * sin(a) + z * cos(a)]
end

# The implementation of the 3-node triangle shell element is in this module. We
# will refer to the functions that we need from this module by referencing them
# relative to the module name.
t3ffm = FEMMShellT3FFModule

# The shell elements have a type `FESetShellT3`. The type of the mesh is the
# plain isoparametric triangle, which is instructed to delegate to the shell
# element.
sfes = FESetShellT3()
accepttodelegate(fes, sfes)

# Use the convenience function to make an instance of the finite element model
# machine for the T3FF shell.
femm = t3ffm.make(IntegDomain(fes, TriRule(1), params.t), mater)

# Construct the requisite fields, geometry and displacement.
# Initialize configuration variables. Displacements are all zero,
# the rotation matrices are all identities.
geom0 = NodalField(fens.xyz)
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
Rfield0 = initial_Rfield(fens)
dchi = NodalField(zeros(size(fens.xyz, 1), 6))

# Apply supports. The clamped end is selected based on the X coordinate.
# We will select all nodes within a box, and the box is slightly inflated using
# the geometrical tolerance based on the spacing of the nodes.

tolerance = min(W / nW, L / nL) / 100

l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
for i in 1:6
    setebc!(dchi, l1, true, i)
end
applyebc!(dchi)
numberdofs!(dchi);

# Associate the finite element model machine with geometry. The shell
# formulation requires knowledge of the normals to the shell surface, and that
# information is computed from the positions of the nodes.

t3ffm.associategeometry!(femm, geom0)

# Assemble the system stiffness matrix.
K = t3ffm.stiffness(femm, geom0, u0, Rfield0, dchi);

# Though load is a concentrated force applied at the center of the beam. First
# we select the node.
nl = selectnode(fens; box = Float64[L L 0 0 0 0], tolerance = tolerance)
# Next we create a mesh of the loaded boundary (i.e. the selected node), so that
# we can integrate along it.
loadbdry = FESetP1(reshape(nl, 1, 1))
# Now we create a finite element model machine that knows how to integrate.
# Since our boundary consists of a single point, we use a point rule.
lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
# Now we create a force intensity to represent the loading.
v = FFlt[0, 0, 0, 0, 0, 0]
v[params.dir] = params.force
fi = ForceIntensity(v);
# Finally we computed the load vector corresponding to the force intensity.
F = distribloads(lfemm, geom0, dchi, fi, 3);

# The system of linear algebraic equations of balance is solved for the displacements and rotations.
U = K \ F
# The vector of unknowns is now distributed into a field.
scattersysvec!(dchi, U[:])

# The deflection in the correct direction at the loaded node is now extracted.
tipdefl = dchi.values[nl, params.dir][1]
@show tipdefl / params.uex * 100

# Generate a graphical display of resultants. The resultants only make sense in
# a coordinate system aligned with the shell surface. Here we construct such a
# coordinate system by using the tangent directions to the surface to construct
# the normal, and the vector along the global X axis will serve as the first
# basis vector of the local cartesian coordinate system.
function updatecsys!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    # first the normal
    cross3!(view(csmatout, :, 3), view(tangents, :, 1), view(tangents, :, 2))
    csmatout[:, 3] ./= norm(view(csmatout, :, 3))
    csmatout[:, 1] .= (1.0, 0.0, 0.0)
    cross3!(view(csmatout, :, 2), view(csmatout, :, 3), view(csmatout, :, 1))
    return csmatout
end
ocsys = CSys(3, 3, updatecsys!)
scalars = []
for nc in 1:3
    fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
    # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
    push!(scalars, ("m$nc", fld.values))
end
vtkwrite("twisted_beam-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
scalars = []
for nc in 1:3
    fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
    # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
    push!(scalars, ("n$nc", fld.values))
end
vtkwrite("twisted_beam-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
scalars = []
for nc in 1:2
    fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
    # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
    push!(scalars, ("q$nc", fld.values))
end
vtkwrite("twisted_beam-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

true


text:
# regular definition
function f(x,y)
  x + y
end

# easier definition
f(x,y) = x + y

g = f

println(g(2, 3))

# convert return type
function l(x, y)::Int8
  return x * y
end;

println(typeof(l(2,3)))

# multiple Return Values
function h(x,y)
  return x+y, x*y
end

println(h(2, 3))

# operators functions
k = +;
println(k(1,2,3,4,5))

# anonymous functions
println(map(x -> x^2 + 2x - 1, [1,3,-1]))

# optional arguments
function m(x, y = 1, z = 3)
  return x + y + z
end
println(m(1))
println(m(1, 3, 5))

# keyword arguments can make these complex interfaces easier to use and extend by
# allowing arguments to be identified by name instead of only by position
# keyword arguments are defined using a semicolon in the signature
function plot(x, y; style="solid", width=1, color="black")
  return x + y + 1, style * color
end
println(plot(1, 2, style = "soft"))
println(plot(2, 3, color = "blue", style = "soft"))

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia