text:
# RANDOMTRANS: function that generates random transactions with dependencies between
#              SKUs based on the values of max_groupsize and max_dependence
function RANDOMTRANS(skus::Int64,
                     orders::Int64,
                     max_group_size::Int64,
                     min_strength::Float64,
                     max_strength::Float64,
                     group_link::Float64,
                     ind_chance::Float64)
    C = Matrix{Float64}(undef,skus,skus) .= 0.0
    D = Matrix{Int64}(undef,skus,skus)   .= 0
    if max_strength > 0.0
        groups = 0
        while groups < skus
            test_correlation = sum(D,dims=1)
            free_skus = findall(==(0),test_correlation)
            group_size = rand(1:max_group_size)
            if group_size <= length(free_skus)
                members = Vector{Int64}(undef,group_size)
                for i = 1:group_size
                    members[i] = rand(1:length(free_skus))
                end
                for i in randperm(skus)[1:rand(1:ceil(Int64,group_size*group_link))]
                    j = rand(1:skus)
                    strength = rand(Uniform(min_strength,max_strength),1)[1]
                    C[i,j] = strength
                    C[j,i] = strength
                    D[i,j] = 1
                    D[j,i] = 1
                end
                for i = 1:group_size
                    for j = 1:group_size
                        if i != j
                            strength = rand(Uniform(min_strength,max_strength),1)[1]
                            local_i = members[i]
                            local_j = members[j]
                            D[local_i,local_j] = 1
                            D[local_j,local_i] = 1
                            C[local_i,local_j] = strength
                            C[local_j,local_i] = strength
                        end
                    end
                end
                groups += group_size
            else
                break
                groups = max_groupmembers
            end
        end
    end
    #transactions = Matrix{Int64}(undef,orders,skus) .= 0
    transactions = spzeros(orders,skus)
    for i = 1:orders
        already_allocated = 0
        skus_order = 1 + floor(abs(rand(Normal(0,3))))
        while already_allocated < skus_order
            new_sku = skus + 1
            while new_sku > skus
                new_sku = 1+ floor(Int64, abs(rand(Normal(0,skus/2))))
            end
            already_allocated += 1
            transactions[i,new_sku] = 1
            if rand() > ind_chance
                for j in randperm(skus)
                    if D[new_sku,j] == 1
                        if transactions[i,j] == 0
                            if already_allocated < skus_order
                                if rand() < C[new_sku,j]
                                    transactions[i,j] = 1
                                    already_allocated += 1
                                end
                            else
                                break
                            end
                        end
                    end
                end
            end
        end
    end
    return transactions
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
using gmsh

function consistent_stl(input::String, output::String; len_scale::Real=1., angle::Real=180.)
    gmsh.initialize()
    gmsh.clear()

    gmsh.merge(input)

    angle = angle * pi / 180
    gmsh.model.mesh.createTopology()
    gmsh.model.mesh.classifySurfaces(angle, true, true)
    gmsh.model.mesh.createGeometry()

    gmsh.model.mesh.set_size(gmsh.model.getEntities(0), len_scale)
    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.optimize("Netgen")

    gmsh.write(output)

    gmsh.clear()
end

ext: jl


text:
export AbstractLearner

using Flux

"""
    (learner::AbstractLearner)(env)

A learner is usually used to estimate state values, state-action values or distributional values based on experiences.
"""
abstract type AbstractLearner end

function (learner::AbstractLearner)(env) end

"""
    get_priority(p::AbstractLearner, experience)
"""
function RLBase.priority(p::AbstractLearner, experience) end

Base.show(io::IO, p::AbstractLearner) =
    AbstractTrees.print_tree(io, StructTree(p), get(io, :max_depth, 10))

function RLBase.update!(
    L::AbstractLearner,
    t::AbstractTrajectory,
    e::AbstractEnv,
    s::AbstractStage,
) end

function RLBase.update!(
    L::AbstractLearner,
    t::AbstractTrajectory,
    e::AbstractEnv,
    s::PreActStage,
)
    update!(L, t)
end


text:
Base.@kwdef mutable struct Options

    # Printing options
    log_verbose::Bool = false
    log_freq::Int = 1000
    timer_verbose::Bool = false
    timer_file::Bool = false
    disable_julia_logger::Bool = true

    # time options
    time_limit::Float64 = 3600_00. #100 hours

    warn_on_limit::Bool = false
    extended_log::Bool = false
    extended_log2::Bool = false
    log_repeat_header::Bool = false

    # Default tolerances
    tol_gap::Float64 = 1e-4
    tol_feasibility::Float64 = 1e-4
    tol_feasibility_dual::Float64 = 1e-4
    tol_primal::Float64 = 1e-4
    tol_dual::Float64 = 1e-4
    tol_psd::Float64 = 1e-7
    tol_soc::Float64 = 1e-7

    check_dual_feas::Bool = false
    check_dual_feas_freq::Int = 1000

    max_obj::Float64 = 1e20
    min_iter_max_obj::Int = 10

    # infeasibility check
    min_iter_time_infeas::Int = 1000
    infeas_gap_tol::Float64 = 1e-4
    infeas_limit_gap_tol::Float64 = 1e-1
    infeas_stable_gap_tol::Float64 = 1e-4
    infeas_feasibility_tol::Float64 = 1e-4
    infeas_stable_feasibility_tol::Float64 = 1e-8

    certificate_search::Bool = true
    certificate_obj_tol::Float64 = 1e-1
    certificate_fail_tol::Float64 = 1e-8

    # Bounds on beta (dual_step / primal_step) [larger bounds may lead to numerical inaccuracy]
    min_beta::Float64 = 1e-5
    max_beta::Float64 = 1e+5
    initial_beta::Float64 = 1.

    # Adaptive primal-dual steps parameters [adapt_decay above .7 may lead to slower convergence]
    initial_adapt_level::Float64 = .9
    adapt_decay::Float64 = .8
    adapt_window::Int64 = 50

    # PDHG parameters
    convergence_window::Int = 200
    convergence_check::Int = 50
    max_iter::Int = 0
    min_iter::Int = 40
    divergence_min_update::Int = 50
    max_iter_lp::Int = 10_000_000
    max_iter_conic::Int = 1_000_000
    max_iter_local::Int = 0 #ignores user setting

    advanced_initialization::Bool = true

    # Linesearch parameters
    line_search_flag::Bool = true
    max_linsearch_steps::Int = 5000
    delta::Float64 = .9999
    initial_theta::Float64 = 1.
    linsearch_decay::Float64 = .75

    # Spectral decomposition parameters
    full_eig_decomp::Bool = false
    max_target_rank_krylov_eigs::Int = 16
    min_size_krylov_eigs::Int = 100
    warm_start_eig::Bool = true
    rank_increment::Int = 1 # 0=multiply, 1 = add
    rank_increment_factor::Int = 1 # 0 multiply, 1 = add

    # eigsolver selection
    #=
        1: Arpack [dsaupd] (tipically non-deterministic)
        2: KrylovKit [eigsolve/Lanczos] (DEFAULT)
    =#
    eigsolver::Int = 2
    eigsolver_min_lanczos::Int = 25
    eigsolver_resid_seed::Int = 1234

    # Arpack
    # note that Arpack is Non-deterministic
    # (https://github.com/mariohsouto/ProxSDP.jl/issues/69)
    arpack_tol::Float64 = 1e-10
    #=
        0: arpack random [usually faster - NON-DETERMINISTIC - slightly]
        1: all ones [???]
        2: julia random uniform (eigsolver_resid_seed) [medium for DETERMINISTIC]
        3: julia normalized random normal (eigsolver_resid_seed) [best for DETERMINISTIC]
    =#
    arpack_resid_init::Int = 3
    arpack_reset_resid::Bool = true # true for determinism
    # larger is more stable to converge and more deterministic
    arpack_max_iter::Int = 10_000
    # see remark for of dsaupd

    # KrylovKit
    krylovkit_reset_resid::Bool = false
    krylovkit_resid_init::Int = 3
    krylovkit_tol::Float64 = 1e-12
    krylovkit_max_iter::Int = 100
    krylovkit_eager::Bool = false
    krylovkit_verbose::Int = 0

    # Reduce rank [warning: heuristics]
    reduce_rank::Bool = false
    rank_slack::Int = 3

    full_eig_freq::Int = 10_000_000
    full_eig_len::Int = 0

    # equilibration parameters
    equilibration::Bool = false
    equilibration_iters::Int = 1000
    equilibration_lb::Float64 = -10.0
    equilibration_ub::Float64 = +10.0
    equilibration_limit::Float64 = 0.9
    equilibration_force::Bool = false

    # spectral norm [using exact norm via svds may result in nondeterministic behavior]
    approx_norm::Bool = true
end


text:
check_list = [
    :h,
    :hbar,
    :ħ,
    :sigmao,
    :sigmax,
    :sigmay,
    :sigmaz,
    :σ₀,
    :σ₁,
    :σ₂,
    :σ₃
]

for each_symbol in check_list
    @test isdefined(each_symbol)
end

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using LinearAlgebra: Symmetric, dot, diagind

@doc raw"""
```julia
restrict_to_inactive!(Q::Symmetric,
                      q::AbstractVector,
                      Z::AbstractVector{BoundConstraint_Enum},
                      lb::AbstractVector{<:Real},
                      ub::AbstractVector{<:Real})
```

```julia
function restrict_to_inactive!(Q::Symmetric,
                               q::AbstractVector,
                               Z::AbstractVector{BoundConstraint_Enum},
                               bc::BoundConstraints{<:Real,1})
```

In-place modification of ``(Q,q)`` that produces
``(\tilde{Q},\tilde{q})`` such that the initial optimization problem:

```math
\tilde{x} = \arg\min \frac{1}{2} x^t Q x + q^tx
```
under these constraints:
```math
x[i] = \left\{\begin{array}{ll}
l[i], & \text{if } Z[i] = -1 \\
u[i], & \text{if } Z[i] = +1
\end{array}\right.
```
is transformed into this linear system:
```math
\tilde{Q}\tilde{x}+\tilde{q}=0
```
"""
function restrict_to_inactive!(Q::Symmetric,
                               q::AbstractVector,
                               Z::AbstractVector{BoundConstraint_Enum},
                               lb::AbstractVector{<:Real},
                               ub::AbstractVector{<:Real})
    Q_ET = eltype(Q)
    n = length(q)

    @assert (n,n)==size(Q)
    @assert n==length(Z)
    @assert n==length(lb)
    @assert n==length(ub)
    @assert Q.uplo=='L' || Q.uplo=='U'

    if Q.uplo=='L'

        for k = 1:n

            if Z[k]!=INACTIVE_BC

                constrained_x_k = ifelse(Z[k]==ACTIVE_LB,
                                         lb[k],
                                         ub[k])

                q[1:(k-1)] .+= constrained_x_k*@view(Q.data[k,1:(k-1)])
                Q.data[k,1:(k-1)] .= zero(Q_ET)

                Q_kk = max(1,abs(Q.data[k,k]))
                q[k] = -constrained_x_k*Q_kk
                Q.data[k,k] = Q_kk

                q[(k+1):n] .+= constrained_x_k*@view(Q.data[(k+1):n,k])
                Q.data[(k+1):n,k] .= zero(Q_ET)

            end # Z[k]!=inactive
        end # k = 1:n
    else
        @assert (Q.uplo=='U') ""

        for k = 1:n

            if Z[k]!=INACTIVE_BC

                constrained_x_k = ifelse(Z[k]==ACTIVE_LB,
                                         lb[k],
                                         ub[k])

                q[1:(k-1)] .+= constrained_x_k*@view(Q.data[1:(k-1),k])
                Q.data[1:(k-1),k] .= zero(Q_ET)

                Q_kk = max(1,abs(Q.data[k,k]))
                q[k] = -constrained_x_k*Q_kk
                Q.data[k,k] = Q_kk

                q[k+1:n] .+= constrained_x_k*@view(Q.data[k,(k+1):n])
                Q.data[k,k+1:n] .= zero(Q_ET)

            end # Z[k]!=inactive
        end # k = 1:n
    end
end
function restrict_to_inactive!(Q::Symmetric,
                               q::AbstractVector,
                               Z::AbstractVector{BoundConstraint_Enum},
                               bc::BoundConstraints{<:Real,1})
    restrict_to_inactive!(Q,q,Z,lower_bound(bc),upper_bound(bc))
end

@doc raw"""
```julia
update_Z!(x::AbstractVector,
          τ::AbstractVector,
          Z::AbstractVector{BoundConstraint_Enum},
          lb::AbstractVector,
          ub::AbstractVector)
```
```julia
update_Z!(x::AbstractVector,
          τ::AbstractVector,
          Z::AbstractVector{BoundConstraint_Enum},
          bc::BoundConstraints{<:Real,1})
```
This function updates `Z` according to `x`, `τ` and bounds `lb`, `ub`
values.

It also count how many changes have be done during this update.

No change means that the algorithm has converged.

**Note:** this function only modifies `Z` and return the number of bad
hypothesis.
"""
function update_Z!(x::AbstractVector,
                   τ::AbstractVector,
                   Z::AbstractVector{BoundConstraint_Enum},
                   lb::AbstractVector,
                   ub::AbstractVector)
    n = length(x)

    @assert n == length(τ)
    @assert n == length(Z)
    @assert n == length(lb)
    @assert n == length(ub)

    count_bad_hypothesis::UInt = 0

    for i in 1:n

        if Z[i]==INACTIVE_BC

            if x[i]<=lb[i]

                count_bad_hypothesis+=1
                Z[i]=ACTIVE_LB

            elseif x[i]>=ub[i]

                count_bad_hypothesis+=1
                Z[i]=ACTIVE_UB

            end

        elseif Z[i]==ACTIVE_LB

            @assert x[i]==lb[i] "Internal error $(x[i]) != $(lb[i])"

            if τ[i]>0
                count_bad_hypothesis+=1
                Z[i]=INACTIVE_BC
            end

        else

            @assert x[i]==ub[i] "Internal error $(x[i]) != $(ub[i])"

            if τ[i]<0
                count_bad_hypothesis+=1
                Z[i]=INACTIVE_BC
            end

        end

    end #  for i in 1:n

    return count_bad_hypothesis
end
function update_Z!(x::AbstractVector,
                   τ::AbstractVector,
                   Z::AbstractVector{BoundConstraint_Enum},
                   bc::BoundConstraints{<:Real,1})
    update_Z!(x,τ,Z,lower_bound(bc),upper_bound(bc))
end

@doc raw"""
```julia
initialize_x_Z(x_init::AbstractArray,
               bc::BoundConstraints)
```

Create (x,Z) from initial guess `x_init` and bound constraints `bc`

`Z` is created by recording how `x_init` fulfils the bound constraints `bc`.

`x` is the projection of `x_init` on the bounded domain [l,b].
"""
function initialize_x_Z(x_init::AbstractArray,
                        bc::BoundConstraints)

    @assert size(x_init) == size(bc)

    x = similar(x_init)
    Z = similar(Array{BoundConstraint_Enum},axes(bc))
    lb = lower_bound(bc)
    ub = upper_bound(bc)

    for (i,(lb_i,x_init_i,ub_i)) in enumerate(zip(lb,x_init,ub))

        if x_init_i<lb_i
            x[i]=lb_i
            Z[i]=ACTIVE_LB
            continue
        end

        if x_init_i>ub_i
            x[i]=ub_i
            Z[i]=ACTIVE_UB
            continue
        end

        x[i]=x_init_i
        Z[i]=INACTIVE_BC
    end

    x, Z
end

@doc raw"""
```julia
update_x!(x::AbstractArray,
          Z::AbstractArray{BoundConstraint_Enum},
          bc::BoundConstraints)
```

Update x value such that:
```math
x[i] = \left\{\begin{array}{ll}
l[i], & \text{if } Z[i] = -1 \\
u[i], & \text{if } Z[i] = +1
\end{array}\right.
```
When ``Z[i]=0`` the ``x[i]`` value is unaffected.
"""
function update_x!(x::AbstractArray,
                   Z::AbstractArray{BoundConstraint_Enum},
                   bc::BoundConstraints)

    @assert size(x) == size(bc) == size(Z)

    for (i,(Z_i,lb_i,ub_i)) in enumerate(zip(Z,lower_bound(bc),upper_bound(bc)))
        if Z_i==ACTIVE_LB
            x[i]=lb_i
            continue
        end

        if Z_i==ACTIVE_UB
            x[i]=ub_i
            continue
        end
    end

    x
end

@doc raw"""
```julia
clean_τ!(τ::AbstractArray{<:Real},
         Z::AbstractArray{BoundConstraint_Enum})
```

By definition τ=-Qx-q. If the algorithm converged, then one must have
τ[i]=0 when the constraint is inactive.

This function updates τ by overwriting τ[i]=0 when Z[i]=inactive.
"""
function clean_τ!(τ::AbstractArray{<:Real},
                  Z::AbstractArray{BoundConstraint_Enum})

    @assert size(τ)==size(Z)

    τ_ELT = eltype(τ)

    for (i,Z_i) in enumerate(Z)
        if Z_i == INACTIVE_BC
            τ[i]=zero(τ_ELT)
        end
    end
    τ
end

# ****************************************************************

# Config stucture
struct Kunisch_Rendl_Conf <: Abstract_BC_QuadSolver_Conf
    _max_iter::Int

    function Kunisch_Rendl_Conf(;
                                max_iter::Int=50)

        @assert max_iter>0

        new(max_iter)
    end
end

max_iter(conf::Kunisch_Rendl_Conf) = conf._max_iter

# ****************************************************************

# Result structure
struct Kunisch_Rendl_Result <: Abstract_BC_QuadSolver_Result
    _cv::Bool
    _iter_count::Int
    _fobj::Real
    _x::AbstractVector{<:Real}
    _τ::AbstractVector{<:Real}

    # For a reason I ignore, using Base.@kwdef leads to an error with
    # Julia 1 By consequence I do not use this macro and write the
    # constructor manually
    function  Kunisch_Rendl_Result(;
                                   _cv,
                                   _iter_count,
                                   _fobj,
                                   _x,
                                   _τ)
        new(_cv,_iter_count,_fobj,_x,_τ)
    end

end
converged(r::Kunisch_Rendl_Result) = r._cv
iteration_count(r::Kunisch_Rendl_Result) = r._iter_count
objective_value(r::Kunisch_Rendl_Result) = r._fobj
multiplier_τ(r::Kunisch_Rendl_Result) = ReadOnlyArray(r._τ)
solution(r::Kunisch_Rendl_Result) = ReadOnlyArray(r._x)

# ****************************************************************

# Kunisch-Rendl algorithm to solve a bound constrained quadratic problem
#
#
# CAVEAT: positive definiteness of Q is no enough to guarantee
# convergence
#
function Kunisch_Rendl(Q::Symmetric{<:Real},
                       q::AbstractVector{<:Real},
                       x_init::AbstractVector{<:Real},
                       bc::BoundConstraints{<:Real,1},
                       conf::Kunisch_Rendl_Conf)
    n = length(q)

    @assert (n,n) == size(Q) "(n,n) == size(Q) $((n,n)) == $(size(Q))"
    @assert n == length(bc)

    (x, Z) = initialize_x_Z(x_init,bc)

    τ = similar(q)
    Q_tilde = similar(Q)
    q_tilde = similar(q)

    has_CV::Bool = false
    local iter_count
    for iter in 1:max_iter(conf)
        iter_count = iter
        Q_tilde .= Q
        q_tilde .= q

        # modify Q,q according to constraints
        restrict_to_inactive!(Q_tilde,q_tilde,Z,bc)

        # x = -Q^{-1}.q
        # TODO: add exception is Q_tilde is too badly conditioned
        try
            x = -Q_tilde\q_tilde
        catch
            @warn "Singular system... Abort..."
            break
        end

        # update x and compute a posteriori multipliers
        update_x!(x,Z,bc)
        τ .= -Q*x-q # for perfs use mul! instead

        # update Z and count bad hypothesis
        count_bad_choice = update_Z!(x,τ,Z,bc)

        # @info "KR: iter=$(_fmt(iter)), wrong_hypothesis=$(_fmt(count_bad_choice))"

        # CV check
        if count_bad_choice==0
            has_CV=true
            break
        end
    end
    # prepare quantities export

    # as τ=-Q*x-q, fobj is 1/2( -τ.x + q.x)
    #
    fobj = (-dot(τ,x) + dot(q,x))/2
    # "true" 0 for inactive constraints
    multiplier_τ = clean_τ!(τ,Z)

    result = Kunisch_Rendl_Result(
        _cv=has_CV,
        _iter_count=iter_count,
        _fobj=fobj,
        _x=x, # no risk as x is local
        _τ=τ # no risk as τ is local
    )

    # @info "Leaving Kunisch-Rendl" result = result

    result
end

# Specialize the generic solve method
#
solve(Q::Symmetric{<:Real},
      q::AbstractVector{<:Real},
      x_init::AbstractVector{<:Real},
      bc::BoundConstraints{<:Real,1},
      conf::Kunisch_Rendl_Conf) = Kunisch_Rendl(Q,q,x_init,bc,conf)


text:
using Base.Test

import FracFin

# write your own tests here
# @test 1 == 2

# grid = 1:10^6

# fBM = FracFin.FractionalBrownianMotion(0.7)
# # sampler = FracFin.CholeskySampler(fBM, grid)
# sampler = FracFin.CircSampler(fBM, grid)
# # sampler = FracFin.CircSampler{FracFin.FractionalBrownianMotion, UnitRange{Int64}}(fBM, grid)
# X = rand(sampler)

# # fGN = FracFin.FractionalGaussianNoise(0.7)
# # sampler = FracFin.CircSampler(fGN, grid)
# # # sampler = CholeskySampler{FractionalGaussianNoise, UnitRange{Int64}}(fGN, grid)
# # X = FracFin.rand(sampler)

# # print(X)

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
using Test
using LyapunovExponents: DEMOS, dimension, solve!, lyapunov_exponents,
    LEProblem, report, objname
using LyapunovExponents.TestTools: test_tangent_dynamics_against_autodiff,
    @test_isapprox_elemwise

@time @testset "Tangent dynamics $(ex.name)" for ex in
        [ex for ex in [f().example for f in DEMOS]
         if ex.tangent_dynamics != nothing]
    num_u = 5  # number of random states at which the systems are tested
    opts = Dict(
        # :verbose => true
    )
    if contains(ex.name, "Linz & Sprott (1999)")
        # linz_sprott_99 is solved with a more accurate and expensive
        # ODE algorithm and with larger t_renorm.  But tangent test
        # (with evolve=true) is done with the default ODE algorithm.
        # So let's old t_renorm as the span.
        # See: [[../src/examples/linz_sprott_99.jl::de_options]]
        opts[:t_evolve] = 1.0
        opts[:t_list] = 0:0.5:1.0
    end
    test_tangent_dynamics_against_autodiff(LEProblem(ex), num_u; opts...)
    if contains(ex.name, "Hénon map")
        # TODO: It fails maybe because some randomly generated states
        # are outside the domain of Hénon map.
        continue
    else
        test_tangent_dynamics_against_autodiff(LEProblem(ex), num_u;
                                               evolve = true,
                                               opts...)
    end
end

@time @testset "Example $(objname(f))" for f in DEMOS
    for dim_lyap in 1:dimension(f().example)
        println()

        demo = f(dim_lyap=dim_lyap)
        @time solve!(demo)
        report(demo)

        sol = demo.solver.sol
        ex = demo.example
        dim = min(dim_lyap, length(ex.known_exponents))

        if contains(ex.name, "standard map") ||
                contains(ex.name, "van der Pol")
            @test_broken sol.converged
        else
            @test sol.converged
        end

        rtol = ex.rtol
        atol = ex.atol
        if contains(ex.name, "standard map")
            rtol *= 20
            atol *= 20
        end

        @test_isapprox_elemwise(
            lyapunov_exponents(sol)[1:dim],
            ex.known_exponents[1:dim],
            rtol = rtol,
            atol = atol,
        )
    end
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# src/affexpr.jl
# Defines all types relating to affine expressions
# - GenericAffExpr              ∑ aᵢ xᵢ  +  c
#   - AffExpr                   Alias for (Float64, Variable)
# - GenericRangeConstraint      l ≤ ∑ aᵢ xᵢ ≤ u
#   - LinearConstraint          Alias for AffExpr
# Operator overloads in src/operators.jl
#############################################################################

#############################################################################
# GenericAffExpr
# ∑ aᵢ xᵢ  +  c
type GenericAffExpr{CoefType,VarType} <: AbstractJuMPScalar
    vars::Vector{VarType}
    coeffs::Vector{CoefType}
    constant::CoefType
end
coeftype{C,V}(::GenericAffExpr{C,V}) = C

Base.zero{C,V}(::Type{GenericAffExpr{C,V}}) = GenericAffExpr{C,V}(V[],C[],zero(C))
Base.one{ C,V}(::Type{GenericAffExpr{C,V}}) = GenericAffExpr{C,V}(V[],C[], one(C))
Base.zero(a::GenericAffExpr) = zero(typeof(a))
Base.one( a::GenericAffExpr) =  one(typeof(a))
Base.copy(a::GenericAffExpr) = GenericAffExpr(copy(a.vars),copy(a.coeffs),copy(a.constant))

# Old iterator protocol - iterates over tuples (aᵢ,xᵢ)
immutable LinearTermIterator{GAE<:GenericAffExpr}
    aff::GAE
end

linearterms(aff::GenericAffExpr) = LinearTermIterator(aff)

Base.start(lti::LinearTermIterator) = 1
Base.done( lti::LinearTermIterator, state::Int) = state > length(lti.aff.vars)
Base.next( lti::LinearTermIterator, state::Int) = ((lti.aff.coeffs[state], lti.aff.vars[state]), state+1)

# More efficient ways to grow an affine expression
# Add a single term to an affine expression
function Base.push!{C,V}(aff::GenericAffExpr{C,V}, new_coeff::C, new_var::V)
    push!(aff.coeffs, new_coeff)
    push!(aff.vars, new_var)
    aff
end
# Add an affine expression to an existing affine expression
function Base.append!{C,V}(aff::GenericAffExpr{C,V}, other::GenericAffExpr{C,V})
    append!(aff.vars, other.vars)
    append!(aff.coeffs, other.coeffs)
    aff.constant += other.constant
    aff
end
# For consistency, allow appending constants and individual variables
Base.append!{C}(aff::GenericAffExpr{C,C}, other::C) = error() # for ambiguity
function Base.append!{C,V}(aff::GenericAffExpr{C,V}, other::C)
    aff.constant += other
    aff
end
Base.append!{C,V}(aff::GenericAffExpr{C,V}, other::V) = push!(aff,one(C),other)

function Base.isequal{C,V}(aff::GenericAffExpr{C,V},other::GenericAffExpr{C,V})
    isequal(aff.constant, other.constant)  || return false
    length(aff.vars) == length(other.vars) || return false
    for i in 1:length(aff.vars)
        isequal(aff.vars[i],   other.vars[i])   || return false
        isequal(aff.coeffs[i], other.coeffs[i]) || return false
    end
    return true
end

# Alias for (Float64, Variable), the specific GenericAffExpr used by JuMP
typealias AffExpr GenericAffExpr{Float64,Variable}
AffExpr() = zero(AffExpr)

Base.isempty(a::AffExpr) = (length(a.vars) == 0 && a.constant == 0.)
Base.convert(::Type{AffExpr}, v::Variable) = AffExpr([v], [1.], 0.)
Base.convert(::Type{AffExpr}, v::Real) = AffExpr(Variable[], Float64[], v)

# Check all coefficients are finite, i.e. not NaN, not Inf, not -Inf
function assert_isfinite(a::AffExpr)
    coeffs = a.coeffs
    for i in 1:length(a.vars)
        isfinite(coeffs[i]) || error("Invalid coefficient $(coeffs[i]) on variable $(a.vars[i])")
    end
end

setobjective(m::Model, sense::Symbol, x::Variable) = setobjective(m, sense, convert(AffExpr,x))
function setobjective(m::Model, sense::Symbol, a::AffExpr)
    if isa(m.internalModel, MathProgBase.AbstractNonlinearModel)
        # Give the correct answer when changing objectives in an NLP.
        # A better approach would be to update and reuse the evaluator
        m.internalModelLoaded = false
    end
    setobjectivesense(m, sense)
    m.obj = convert(QuadExpr,a)
end

# Copy an affine expression to a new model by converting all the
# variables to the new model's variables
function Base.copy(a::AffExpr, new_model::Model)
    AffExpr(copy(a.vars, new_model), copy(a.coeffs), a.constant)
end

function getvalue(a::AffExpr)
    ret = a.constant
    for it in 1:length(a.vars)
        ret += a.coeffs[it] * getvalue(a.vars[it])
    end
    ret
end

##########################################################################
# GenericRangeConstraint
# l ≤ ∑ aᵢ xᵢ ≤ u
# The constant part of the internal expression is assumed to be zero
type GenericRangeConstraint{TermsType} <: AbstractConstraint
    terms::TermsType
    lb::Float64
    ub::Float64
end

#  b ≤ expr ≤ b   →   ==
# -∞ ≤ expr ≤ u   →   <=
#  l ≤ expr ≤ ∞   →   >=
#  l ≤ expr ≤ u   →   range
function sense(c::GenericRangeConstraint)
    if c.lb != -Inf
        if c.ub != Inf
            if c.ub == c.lb
                return :(==)
            else
                return :range
            end
        else
                return :(>=)
        end
    else #if c.lb == -Inf
        c.ub == Inf && error("'Free' constraint sense not supported")
        return :(<=)
    end
end

function rhs(c::GenericRangeConstraint)
    s = sense(c)
    s == :range && error("Range constraints do not have a well-defined RHS")
    s == :(<=) ? c.ub : c.lb
end

# Alias for AffExpr
typealias LinearConstraint GenericRangeConstraint{AffExpr}

function Base.copy(c::LinearConstraint, new_model::Model)
    return LinearConstraint(copy(c.terms, new_model), c.lb, c.ub)
end

function addconstraint(m::Model, c::LinearConstraint)
    push!(m.linconstr,c)
    if m.internalModelLoaded
        if method_exists(MathProgBase.addconstr!, (typeof(m.internalModel),Vector{Int},Vector{Float64},Float64,Float64))
            assert_isfinite(c.terms)
            indices, coeffs = merge_duplicates(Cint, c.terms, m.indexedVector, m)
            MathProgBase.addconstr!(m.internalModel,indices,coeffs,c.lb,c.ub)
        else
            Base.warn_once("Solver does not appear to support adding constraints to an existing model. Hot-start is disabled.")
            m.internalModelLoaded = false
        end
    end
    return LinConstrRef(m,length(m.linconstr))
end
addconstraint(m::Model, c::Array{LinearConstraint}) =
    error("The operators <=, >=, and == can only be used to specify scalar constraints. If you are trying to add a vectorized constraint, use the element-wise dot comparison operators (.<=, .>=, or .==) instead")

function addVectorizedConstraint(m::Model, v::Array{LinearConstraint})
    ret = Array(LinConstrRef, size(v))
    for I in eachindex(v)
        ret[I] = addconstraint(m, v[I])
    end
    ret
end


text:
# An infinite number of mathematicians walk into a bar...
using Colors
using Compat

colors = distinguishable_colors(9)

box(w, h, n) =
    empty |> fillcolor(colors[n % 9 + 1]) |> size(w, h)

cut(w, h, ::@compat Tuple{Escher.Vertical, Bool}) =
    (w, h/2)
cut(w, h, ::Any) =
    (w/2, h)

stack(w, h, d, n) =
    n == 0 ?
        empty :
        flow(d[1][1], [ box(cut(w, h, d[1])..., n),
            stack(cut(w, h, d[1])..., vcat(d[2:end], d[1]), n-1)], reverse=d[1][2])

directions = [
    (vertical, false),
    (horizontal, true),
    (vertical, true),
    (horizontal, false)
]

main(win) = inset(Escher.middle,
             size(100vw, 100vh, empty),
             stack(80vmin, 80vmin, directions, 24) |> paper(2))


text: ImportToolbox(str) = arcpy.ImportToolbox(str)


text:
using ApproxFun, SingularIntegralEquations, RiemannHilbert, Test
using RiemannHilbert.KdV

V = x -> 1.2exp(-x^2)
R = ReflectionCoefficient(V)
@test R(0.1) ≈ -0.9651837057639975 + 0.06801873361192057im  # From ISTPackage
@test R(0.01) ≈ -0.9996431456648749 + 0.0070094571751774535im
@test R(0.0) ≈ -1.0
@test R(-0.1) ≈ conj(R(0.1))

k = 1im
a,x₀ = endpoints(domain(R.V₋))
b = rightendpoint(domain(R.V₊))
D = Derivative()
V₋,V₊ = R.V₋, R.V₊

kk = range(-1,1,length=10) .+ range(0,1,length=10)'*im

ψ = (D,a,x,k) -> begin
    ψ = [ivp(); D^2  + (V₋ + k^2)] \ [exp(-im*k*a), -im*k*(exp(-im*k*a)), 0.0]
    ψ(x)
end

kk = range(-10,10,length=300)' .+ range(0,5,length=300)*im
f = ψ.(Ref(D),a,0.0,kk)

portrait(f, PTstepmod)

ψ(D,a,0.0,10im)
using ComplexPhasePortrait
portrait(f)

@time ρ = tFun(R, -5.0..5, 300)
@test ρ(0.1) ≈ R(0.1)
k = Fun(identity, space(ρ))
G = (t,x) -> [1-abs2.(ρ) -conj.(ρ)*exp(-2im*k*x-8im*k^3*t);
        ρ*exp(2im*k*x+8im*k^3*t)           1.0]

Gx = (t,x) -> [0 2im*k*conj.(ρ)*exp(-2im*k*x-8im*k^3*t);
        2im*k*ρ*exp(2im*k*x+8im*k^3*t)           0.0]

# Φ⁺ = Φ⁻*G
# Φₓ⁺ - Φₓ⁻*G = Φ⁻*Gₓ
# Φₓ = ψ*Φ
# ψ₊ - ψ₋ = Φ⁻*Gₓ*inv(Φ₊)

# U - (G-I)*C₋U = G-I
#

import RiemannHilbert: fpstieltjesmatrix
import ApproxFun: transform
@time Φ = transpose(rhsolve(transpose(G), 2*4*200))
h = 0.0001
@time Φ_h = transpose(rhsolve(transpose(G(0.0,h)), 2*4*200))
n = 2*4*200; S₋ = fpstieltjesmatrix(space(U)[1,1], n÷2, n÷2)
(z = 0.1+eps()im; (Φ_h(z)-Φ(z))/h) - (z = 0.1-eps()im; ((Φ_h(z)-Φ(z))/h)*G(0.,0.)(0.1))

Φ_h(0.1+eps()im)-Φ_h(0.1-eps()im)*G(0.0,h)(0.1)

(Φ₋*Gx(0,0))(0.1)

cauchy(Φ₋*Gx(0,0)*quickinv(Φ₊), z)*Φ(z)

(z=0.1+eps()im; cauchy(Φ₋*Gx(0,0)*quickinv(Φ₊), z)*Φ(z)) -
    (z=0.1-eps()im; cauchy(Φ₋*Gx(0,0)*quickinv(Φ₊), z)*Φ(z))*G(0,0)(0.1) -
            (Φ₋*Gx(0,0))(0.1)

(istieltjes(Φ)/(-2π*im)+Φ₋)(0.1)

Φ₊ = (istieltjes(Φ)*(-2π*im)+Φ₋)
Φ(0.1-eps()im)

Φ₊(0.1) - Φ₋(0.1)

Φ(0.1+eps()*im) - Φ(0.1-eps()*im)

I+stieltjes(istieltjes(Φ),0.1+eps()*im) - Φ(0.1+eps()*im)

I+stieltjes(istieltjes(Φ),0.1-eps()*im) - Φ(0.1-eps()*im)
stieltjes(istieltjes(Φ),0.1+eps()*im)-stieltjes(istieltjes(Φ),0.1-eps()*im)

istieltjes(Φ)(0.1)*(-2π*im)

istieltjes(Φ)(0.1)

Φ₋(0.1) - Φ(0.1-eps()*im)

U(0.1)/(-2π*im)

Fun(I+[Fun(space(U)[k,j], ApproxFun.transform(space(U)[k,j], S₋*U[k,j].coefficients)) for k=1:2,j=1:2]) |> space

Φ₋(0.1)
Φ(0.1-eps()im)

V(0.1

Φ(0.1+eps()im) - Φ(0.1-eps()im)*G(0,0)(0.1)

istieltjes(Φ)(0.1)*(-2π*im) - (Φ(0.1+eps()im) - Φ(0.1-eps()im))

G⁻¹(0.1) - inv(G(0,0)(0.1))

G⁻¹(0.1)*Φ(0.1-eps()im) - Φ(0.1+eps()im)

U = istieltjes(Φ)
Φ₋ =  Fun(Fun(I+[Fun(space(U)[k,j], ApproxFun.transform(space(U)[k,j], S₋*U[k,j].coefficients)) for k=1:2,j=1:2]),
            ApproxFun.ArraySpace(Chebyshev(domain(ρ)), 2, 2))

domain(ρ)
Φ₋⁻¹ = quickinv(Φ₋)
G⁻¹ = quickinv(G)

z = 190000.0im; 2im*(z*[1,1]'cauchy(Φ₋*Gx(0,0)*quickinv(Φ₊), z)*Φ(z))[1]

cauchy(

V(0.0)
sech(0.0)

space(Gx)

inv.(G)

G[1,1]*Φ₋⁻¹[1,1]

ncoefficients(G[1,1])
ncoefficients(Φ₋⁻¹)

space(Φ₋⁻¹)

M = [Fun(sp,transform(sp,[Vi[p][k,j] for p=1:length(Vi)])) for k=1:2, j=1:2]
k=j=1;Fun(sp,transform(sp,[Vi[p][k,j] for p=1:length(Vi)])) |> typeof
typeof(sp)
ϕ₋⁻¹(0.1) - inv(Φ(0.1-eps()im))

S₋*U[1,1].coefficients

inv.(Φ₋)

@profiler stieltjes(U,0.1+eps()im)

Φ(0.1-eps()*im)

U = istieltjes(Φ)
@time Φ(0.1⁺)
@profiler stieltjes(U,0.1+0.1im)

Fun(x -> Φ((x)⁺), domain(ρ), 20)

ncoefficients(ρ)

(I+stieltjes(U,1+im)) - Φ(1+im)

([1 1]*Φ(100.0im))[1]

norm(Φ(0.1+eps()im) - Φ(0.1-eps()im)*G(0.1))

hilbert(

D^2 + Fun(V, -10..10)

z = Fun(ℂ)
-2im*(z*Φ[2,1])(Inf)
-2im*1000Φ[2,1](1000)
0.1sech(0.0)

Φ(100.0im)

Φ(0.1-eps()im)*G(0.1)

domain(ρ)

v = tvalues(R, Chebyshev(-10..10), 100)

R(8.0)

plot(real.(v))
    plot!(imag.(v))

R(0.0)
R(-0.0000)

ρ(0.1)
R(0.1)

plot(ρ)

ρ.coefficients

ρ.coefficients

f = R
n = 100
d = Chebyshev(-10..10)

@time R(1.0)
ret = Vector{ComplexF64}(undef,100)
@time Threads.@threads for k=1:length(ret)
    ret[k] = R(k/10)
end

@time for k=1:100
    ret[k] = R(k/10)
end

@time pmap(R, 1.0:4)

ρ = Fun(R, -10..10, 200)

@time R(10.0)

ks = -10.0:0.1:10.0
Rks = R.(ks)

plot(ks, real.(Rks))
    plot!(ks, imag.(Rks))

R(0.1+0.1im)

#####
# Factor KdV
######


text:
using Test
using P4est

@testset "P4est.jl tests" begin
  @testset "p4est_connectivity_new_periodic" begin
    @test p4est_connectivity_new_periodic() isa Ptr{p4est_connectivity}
  end

  connectivity = p4est_connectivity_new_periodic()
  @testset "p4est_connectivity_is_valid" begin
    @test p4est_connectivity_is_valid(connectivity) == 1
  end

  @testset "unsafe_wrap" begin
    obj = unsafe_wrap(connectivity)
    @test obj.num_vertices == 4
  end
end


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

struct SubnetsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a subnet in the specified virtual network.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: subnetParameters::Subnet (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Subnet
"""
function subnetsCreateOrUpdate(_api::SubnetsApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, subnetParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", Subnet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}", ["azure_auth"], subnetParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified subnet.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function subnetsDelete(_api::SubnetsApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified subnet by virtual network and resource group.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: Subnet
"""
function subnetsGet(_api::SubnetsApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Subnet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all subnets in a virtual network.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SubnetListResult
"""
function subnetsList(_api::SubnetsApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SubnetListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export subnetsCreateOrUpdate, subnetsDelete, subnetsGet, subnetsList


text:
abstract type RandomSearch end
include("purerandomsearch.jl")
include("simulatedannealing.jl")
include("particleswarm.jl")


text:
#=
This is a multiline comment
with two lines =#

using PyPlot

x = Int[]
for i in -3:3
    push!(x, i)
end

y1 = x.^2
y2 = x.^1

scatter(x, y1)
plot(x, y2)

filename = tempname()
savefig(filename)

fig2 = figure()
ax = subplot(111)

ax.plot(x, -y1)
ax.scatter(x, 5 .* rand(length(x)))

fig2.savefig(filename)


text:
"""
    BAEstep!(lnL, ind, nod, w, U, pars, Nx, Q)
Compute E step proposed by Bock and Aitkin (1981). The integral is apploximated by discrete normal quadrature.
This function id in-place type.

Non adaptive, non GH version

# Arguments

`lnL` A p(n of subjects) × l(n of nodes) matrix of log likelohood.
`post` A p × l matrix of the apprroximated posterior density.
`pind` Array of arrays that contain the locations of person in each groupes.
`iind` Array of arrays that contains the item location of each records.
`gh` Vector of DNq.
`U` A matrix of item responses.
`pars` Array of IRT model parameters.
`Q` Q matrix.
"""
function BAEstep!(lnL, post, pind, iind, gh::Array{DNq}, U, pars, Q)
    # Reset
    for i in axes(lnL, 1), j in axes(lnL, 2)
            @fastmath @inbounds lnL[i, j] = 0.0
    end
    for g in axes(pind, 1)
        @fastmath @inbounds Threads.@threads for p in pind[g] # person index (Vector of Vector)
            for i in iind[p] # item index (Vector of Vector)
                if U[p,i] !== missing
                    for l::Int64 in axes(gh[g].w, 1)
                        lnL[p,l] += @views log(irf(pars[i], gh[g].n[l,:], U[p,i], Q[i,:]))
                    end
                end # of l
            end # of item
            # approximate posterior distribution
            P = @fastmath @views exp.(lnL[p,:]) .* gh[g].w
            post[p,:] = @fastmath P ./ sum(P)
        end # of person
    end # of group
    return nothing
end

function initializematrix(P, I; use_static = false)
    if use_static
        return @MMatrix zeros(Float64, P, I)
    else
        return zeros(Float64, P, I)
    end
end


text:
const VarCostArgs = Union{Float64, NTuple{2, Float64}, Vector{NTuple{2, Float64}}}

abstract type OperationalCost <: DeviceParameter end

mutable struct VariableCost{T}
    cost::T
end

get_cost(vc::PowerSystems.VariableCost) = vc.cost
Base.length(vc::PowerSystems.VariableCost) = length(vc.cost)
Base.getindex(vc::PowerSystems.VariableCost, ix::Int64) = getindex(vc.cost, ix)


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

using PyCall
p = pyimport("pyplasm")

geom_0 = hcat([[x] for x=0.:10]...)
topol_0 = [[i,i+1] for i=1:9]
geom_1 = hcat([[0.],[1.],[2.]]...)
topol_1 = [[1,2],[2,3]]

model_0 = (geom_0,topol_0)
model_1 = (geom_1,topol_1)
model_2 = Lar.larModelProduct(model_0, model_1)
model_3 = Lar.larModelProduct(model_2, model_1)

p["VIEW"](Plasm.lar2hpc(model_3...))
V,CV = model_3

hpc = Plasm.hpc_exploded( (V,[CV]) )(1.5,1.5,1.5)
Plasm.view(hpc)

hpc = Plasm.lar2exploded_hpc(model_3[1], model_3[2])(1.5,1.5,1.5)
Plasm.view(hpc)

shape = [10,10,2]
cubes = Lar.larCuboids(shape,true)
V,FV = cubes[1],cubes[2][3]
Plasm.viewexploded(V,FV)()

V,EV = Lar.larCuboids([30])
V1 = hcat( map(u->[u; 0.0], V )...)
W = hcat( map(u->[cos(u*(2π/30)); sin(u*(2π/30))], V )...)
Plasm.viewexploded(V1,EV)()
Plasm.viewexploded(W,EV)()

V,FV = Lar.larCuboids([30,6])
V2 = [2π/30 0; 0 1/6] * V
W = [V2[:,k] for k=1:size(V,2)]
Z = hcat( map(p->let (u,v) = p; [v*cos(u); v*sin(u)] end, W) ...)
Plasm.viewexploded(V2,FV)()
Plasm.viewexploded(Z,FV)()

V,CV = Lar.larCuboids([16,32,1])
V3 = [π/16 0 0; 0 2π/32 0; 0 0 1] * V
W = [V3[:,k]-[π/2,π,1/2] for k=1:size(V3,2)]
Z = hcat( map(p->let (u,v,w) = p; [w*cos(u)*cos(v); w*cos(u)*sin(v); w*sin(u)] end, W) ...)
Plasm.viewexploded(V3,CV)()
Plasm.viewexploded(Z,CV)(2,2,2)


text:
#RSDeltaSigma: Quantizer functionnality
#-------------------------------------------------------------------------------

#Version from simulateDSM.m
#-------------------------------------------------------------------------------
#VERIFYME: Did this one actually get called? Was the other version the one being called??
"""`ds_quantize_DSM(y,n)`

Quantize y to:
 - an odd integer in [-n+1, n-1], if n is even, or
 - an even integer in [-n, n], if n is odd.

This definition gives the same step height for both mid-rise
and mid-tread quantizers.
"""
function ds_quantize_DSM(y, n::Int)
	if mod(n,2)==0 #mid-rise quantizer
		v = 2*floor.(0.5*y) .+ 1
	else #mid-tread quantizer
		v = 2*floor.(0.5*(y .+ 1))
	end

	#Limit the output
	for qi in 1:length(n) #Loop for multiple quantizers
		L = n[qi]-1
		i = v[qi,:] .> L
		if any(i)
			v[qi,i] .= L
		end
		i = v[qi,:] .< -L
		if any(i)
			v[qi,i] .= -L
		end
	end
	return v
end
#Pretty sure mod(n,2) code needs to be adapted for Vector:
ds_quantize_DSM(y, n::Vector{Int}) =
	throw(ErrorException("ds_quantizer: does not support n::Vector"))

#Version from ds_quantize.m:
#-------------------------------------------------------------------------------
#=
function v = ds_quantize(y,n)
%v = ds_quantize(y,n=2)
%Quantize y to
% an odd integer in [-n+1, n-1], if n is even, or
% an even integer in [-n+1, n-1], if n is odd.
%This definition gives the same step height for both mid-rise
%and mid-tread quantizers.
%n can be a column vector which specifies how to quantize the rows of y
if nargin<2
    n=2;
end
if length(n)==1
    n = n*ones(size(y));
elseif size(n,1)==size(y,1)
    n = n*ones(1,size(y,2));
end

i = rem(n,2)==0;
v = zeros(size(y));
v(i) = 2*floor(0.5*y(i))+1;	% mid-rise quantizer
v(~i) = 2*floor(0.5*(y(~i)+1));	% mid-tread quantizer

% Limit the output
L = n-1;
for m=[-1 1]
    i = m*v>L;
    if any(i(:))
	v(i) = m*L(i);
    end
end
=#

#Last line


text:
module Utils

using DataFrames
using CSV
using Dates
using TSML
import AutoMLPipeline
AMLP = AutoMLPipeline

export readdata,groupedmean,extractgrmean,hourly,extractsensor,FishDType,FVal

# create dynamic dispatch based on data type: abm, cad2, cad4
abstract type FishDType end

struct FVal{x} <: FishDType
end

FVal(x) = FVal{x}()

function groupedmean(tpe::FishDType,dv::Union{DataFrameRow,Vector}) end
function extractgrmean(tpe::FishDType,data::DataFrame)  end
# end of dynamic dispatch

# extract temp/salinity/oxygen
function extractsensor(fname,sy::Symbol)
  tempdata=readdata(fname,"y-m-d H:M:S+S:S",:observed_timestamp)
  if :label in names(tempdata)
	 select!(tempdata,Not(:label)) # remove redundant column
  end
  if ncol(tempdata) == 3
     hourlytemp=hourly(tempdata[:,2:3])
  else
     hourlytemp=hourly(tempdata)
  end
  rename!(hourlytemp,Dict(:Value => sy))
  return hourlytemp
end

function readdata(file::String,datefmt::String="y-m-d H:M:S",datefield::Symbol=:Timestamp)
	dataloc = joinpath(@__DIR__,"../data/"*file)
	data = CSV.File(dataloc) |> DataFrame
	data[!,datefield] = DateTime.(data[!,datefield],datefmt)
	return data
end

function groupedmean(::FVal{:abm},dv::Vector)
	#row-wise impute
	v = interp(dv) |> locf() |> nocb()
	len = length(v)
	wndx = collect(0:0.5:len)[2:2:2len] # depths midpoints
	@assert length(v) == length(wndx)
	return sum(wndx .* v)/sum(v)
end

function groupedmean(::FVal{:cad2},dv::Vector)
	#row-wise impute
	v = interp(dv) |> locf() |> nocb()
	wndx = collect(0.1:0.2:10.5) # depths midpoints
	@assert length(v) == length(wndx)
	return sum(wndx .* v)/sum(v)
end

function groupedmean(::FVal{:cad4},dv::Vector)
	#row-wise impute
	v = interp(dv) |> locf() |> nocb()
	wndx = collect(0.1:0.2:8.7) # depths midpoints
	@assert length(v) == length(wndx)
	return sum(wndx .* v)/sum(v)
end

function groupedmean(::FVal{:scot},dv::Vector)
	#row-wise impute
	v = interp(dv) |> locf() |> nocb()
	wndx = collect(0.1:0.2:2.7) # depths midpoints
	@assert length(v) == length(wndx)
	return sum(wndx .* v)/sum(v)
end

function groupedmean(tpe::FishDType,dv::DataFrameRow)
	groupedmean(tpe,Vector(dv))
end

function groupedmean(tpe::FishDType,data::DataFrame)
	[groupedmean(tpe,x) for x in eachrow(data)]
end

function extractgrmean(tpe::FVal{:abm},data::DataFrame)
	# check specific data
	@assert eltype(data[:,1]) <: DateTime
	@assert ncol(data) == 23
	# global impute
	vals = float.(data[:,2:(end-2)]) # extract depths
	# columnwise impute
	vals = interp(vals) |> locf() |> nocb()
	fvals = groupedmean(tpe,vals)
	cdata = DataFrame(time=data.Timestamp,mean=fvals)
end

function extractgrmean(tpe::Union{FVal{:cad2},FVal{:cad4},FVal{:scot}},data::DataFrame)
	# check specific data
	@assert eltype(data[:,1]) <: DateTime
	#@assert ncol(data) == 54
	# global impute
	vals = float.(data[:,2:end]) # extract depths
	# columnwise impute
	vals = interp(vals) |> locf() |> nocb()
	fvals = groupedmean(tpe,vals)
	cdata = DataFrame(time=data.Timestamp,mean=fvals)
end

function hourly(df::DataFrame,plot::Bool=false)
	valnner = DateValNNer(Dict(:strict=>false,:nnsize=>2))
	valgator = DateValgator()
	valizer = DateValizer()
	outliernicer = Outliernicer()
	pltr = Plotter()
	if plot == true
		pipeline = @pipeline valgator |> valnner |> pltr
	else
		pipeline = @pipeline valgator |> valnner
	end
	fit!(pipeline,df)
	ct=TSML.transform!(pipeline,df)
	ct
end

end

max_stars_repo_name: IBM/precision-aquaculture

max_issues_repo_name: IBM/precision-aquaculture

max_forks_repo_name: IBM/precision-aquaculture


text:
"""
            Types needed for gadget read routines.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

mutable struct Header
    npart::Vector{Int32}                # an array of particle numbers per type in this snapshot
    massarr::Vector{Float64}            # an array of particle masses per type in this snapshot - if zero: MASS block present
    time::Float64                       # time / scale factor of the simulation
    z::Float64                          # redshift of the simulation
    flag_sfr::Int32                     # 1 if simulation was run with star formation, else 0
    flag_feedback::Int32                # 1 if simulation was run with stellar feedback, else 0
    nall::Vector{UInt32}                # total number of particles in the simulation
    flag_cooling::Int32                 # 1 if simulation was run with cooling, else 0
    num_files::Int32                    # number of snapshots over which the simulation is distributed
    boxsize::Float64                    # total size of the simulation box
    omega_0::Float64                    # Omega matter
    omega_l::Float64                    # Omega dark enery
    h0::Float64                         # little h
    flag_stellarage::Int32              # 1 if simulation was run with stellar age, else 0
    flag_metals::Int32                  # 1 if simulation was run with metals, else 0
    npartTotalHighWord::Vector{UInt32}  # weird
    flag_entropy_instead_u::Int32       # 1 if snapshot U field contains entropy instead of internal energy, else 0
    flag_doubleprecision::Int32         # 1 if snapshot is in double precision, else 0
    flag_ic_info::Int32
    lpt_scalingfactor::Float32
    fill::Vector{Int32}                 # the HEAD block needs to be filled with zeros to have a size of 256 bytes

    function Header(npart::Vector{Int32}=Int32.([0,0,0,0,0,0]),
           massarr::Vector{Float64}=zeros(6),
           time::Float64=0.,
           z::Float64=0.,
           flag_sfr::Int32=Int32(0),
           flag_feedback::Int32=Int32(0),
           nall::Vector{UInt32}=UInt32.([0,0,0,0,0,0]),
           flag_cooling::Int32=Int32(0),
           num_files::Int32=Int32(0),
           boxsize::Float64=0.,
           omega_0::Float64=0.,
           omega_l::Float64=0.,
           h0::Float64=0.,
           flag_stellarage::Int32=Int32(0),
           flag_metals::Int32=Int32(0),
           npartTotalHighWord::Vector{UInt32}=UInt32.([0,0,0,0,0,0]),
           flag_entropy_instead_u::Int32=Int32(0),
           flag_doubleprecision::Int32=Int32(0),
           flag_ic_info::Int32=Int32(0),
           lpt_scalingfactor::Float32=Float32(0.),
           fill::Vector{Int32}=Int32.(zeros(12)))

          new(npart,
              massarr,
              time,
              z,
              flag_sfr,
              flag_feedback,
              nall,
              flag_cooling,
              num_files,
              boxsize,
              omega_0,
              omega_l,
              h0,
              flag_stellarage,
              flag_metals,
              npartTotalHighWord,
              flag_entropy_instead_u,
              flag_doubleprecision,
              flag_ic_info,
              lpt_scalingfactor,
              fill)
    end
end

mutable struct Info_Line
    block_name::String              # name of the data block, e.g. "POS"
    data_type::DataType             # datatype of the block, e.g. Float32 for single precision, Float64 for double
    n_dim::Int32                    # number of dimensions of the block, usually 1 or 3
    is_present::Vector{Int32}       # array of flags for which particle type this block is present,
                                    # e.g. gas only:  [ 1, 0, 0, 0, 0, 0 ]
                                    # e.g. gas + BHs: [ 1, 0, 0, 0, 0, 1 ]

    function Info_Line(block_name="", data_type=Float32, n_dim=Int32(0),
                        is_present=Int32.(zeros(6)))

        new(block_name, data_type, n_dim, is_present)
    end
end

mutable struct part_single
      ptype::Int64
      pos::Vector{Float32}
      vel::Vector{Float32}
      id::UInt32
      m::Float32
      U::Float32
      rho::Float32
      sml::Float32

      function part_single(ptype=-1, pos=Float32.(zeros(3)), vel=Float32.(zeros(3)),
                        id=UInt32(0), m=Float32(0.), U=Float32(0.), rho=Float32(0.),
                        sml=Float32(0.))

             new(ptype, pos, vel, id, m, U, rho, sml)

      end
end

mutable struct part_double
      ptype::Int64
      pos::Vector{Float64}
      vel::Vector{Float64}
      id::UInt32
      m::Float64
      U::Float64
      rho::Float64
      sml::Float64

      function part_double(ptype=-1, pos=zeros(3), vel=zeros(3),
                        id=UInt32(0), m=0., U=0., rho=0.,
                        sml=0.)

             new(ptype, pos, vel, id, m, U, rho, sml)

      end
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


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

using Markdown
using InteractiveUtils

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

# ╔═╡ a7ef92d7-41c1-4f7d-a0f9-b8037237af4f
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add(url="https://github.com/mkitti/GaussianBeamLens.git")
	Pkg.add("GR")
	Pkg.add("PlutoUI")
end

# ╔═╡ bca271a4-ee26-48ab-bef8-e6c28ae76351
begin
	using GaussianBeamLens
	using GR
	using PlutoUI
end

# ╔═╡ 10e09fe7-cc23-4487-b1dc-b0a545bf2489
begin
	#cm = GR.COLORMAP_INFERNO
	colormaps = [
			string(getfield(GR,name)) => replace(string(name),"COLORMAP_" => "")
			for name in names(GR, all=true)
			if startswith(string(name),"COLORMAP_")
				]
	println(colormaps)
	md"""
	Colormap:
		$(@bind cm Select(colormaps; default = "45"))
	"""
end

# ╔═╡ b1e13be1-d608-43a1-a7d5-29a6c3f98229
md"""
x0:
$(@bind x0 Slider(-240:10:240; default=0, show_value=true))

w0:
$(@bind w0 Slider(2:1:50; default=25, show_value=true))

f:
$(@bind f Slider(0:10:2000; default=1000, show_value=true))

Lens\_z:
$(@bind Lens_z Slider(0:10:2000; default=1000, show_value=true))
"""

# ╔═╡ 4b7b50c6-1304-4c81-b2ac-92ecf7a3d4d1
begin
xlabel("z")
ylabel("x")
GR.heatmap(GaussianBeamLens.z,GaussianBeamLens.x,GaussianBeamLensPropagate(x0, w0, f, Lens_z)'; colormap = parse(Int,cm))
end

# ╔═╡ Cell order:
# ╠═a7ef92d7-41c1-4f7d-a0f9-b8037237af4f
# ╠═bca271a4-ee26-48ab-bef8-e6c28ae76351
# ╠═10e09fe7-cc23-4487-b1dc-b0a545bf2489
# ╠═4b7b50c6-1304-4c81-b2ac-92ecf7a3d4d1
# ╠═b1e13be1-d608-43a1-a7d5-29a6c3f98229

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
#STOCHASTIC MODELING FUNCTIONS FOR MC HAMMER
#by Eric Torkia, August 2020

#Markov Chains, Martingales and other useful stochastic processes

"""
    cmatrix(dimensions::Int64)

cmatrix produces an N x N matrix of 0s where the last row are all 1s.

3×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  1.0  1.0

"""
function cmatrix(dimensions::Int64)

#Create Vector with final value  = 1
    sample_array = zeros(dimensions-1)
    push!(sample_array, 1)

    #Add columns to Matrix
    cmatrix = sample_array

    for i = 1:(dimensions-1)
        cmatrix = hcat(cmatrix, sample_array)
    end
    return cmatrix

end

"""
    markov_a(t_matrix)

markov_a produces the calculated end state of a Markov chain using a sqaure transition matrix.

"""
function markov_a(t_matrix)
    md = size(t_matrix,1)
    IMatrix =  Matrix(I,md,md)
    analytic_markov  =  transpose(t_matrix - IMatrix) + cmatrix(md)
    analytic_markov = analytic_markov^-1
    return      analytic_markov[:,md]
end

#Support function to calculate the Markov chain values at specific point in the chain.
function markov_states(t_matrix, start_arr)
    states = size(t_matrix,1)
    state_mat = t_matrix .* start_arr
    new_state_arr = []

    for state = 1:states
        state_value = sum(state_mat[:,state])
        push!(new_state_arr, state_value)
    end
    return new_state_arr
end

## markov_ts

"""
    markov_ts(t_matrix, start_arr, trials=1)

*markov_ts* is a time-series method that allows you to see the transition states trial by trial.
Extremely useful for market share or account receivables problems where you want to see the states change over time.

*t_matrix* is the transition Matrix
*start_arr* is the starting values for the chain.
 *trials* is the number of iterations you want to run through the Markov Chain Process.

"""
function markov_ts(t_matrix, start_arr, trials=1)
    mts_1 = markov_states(t_matrix, start_arr) #Markov Chaine after 1 trial
    mts_series = mts_1 # setup the first column in the results
    mts_x = mts_1

    for trial = 2:trials
        mts_x = markov_states(t_matrix, mts_x)
        mts_series = hcat(mts_series, mts_x)
    end
    return mts_series
end

## Martingales

"""
    marty(Wager, GamesPlayed; GameWinProb = 0.5, CashInHand = Wager)

*In probability theory, a martingale is a sequence of random variables (i.e., a stochastic process) for which, at a particular time, the conditional expectation
of the next value in the sequence, given all prior values, is equal to the present value.* (Wikipedia)

The marty function is designed to simulate a Martigale such as that everytime a wager is lost, the next bet doubles the wagered amount to negate the previous loss.
The resulting vector is the balance of cash the gambler has in hand at any given point in the Martingale process.

*GameWinProb* is the estimated probability of winning.
*CashInHand* is the starting balance for the martigale. At times, this parameter can make a difference in whether your survive the process or go home broke.

"""
function marty(Wager, GamesPlayed; GameWinProb = 0.5, CashInHand = Wager)
    CurrentWager = Wager
    Balance = CashInHand
    BalResults = []
        for i = 1:GamesPlayed
            win_loose_sim = rand(Bernoulli(GameWinProb))

            if win_loose_sim == true
                    Balance = Balance + CurrentWager
                    CurrentWager = Wager
                elseif win_loose_sim == false
                    Balance = Balance - CurrentWager
                    CurrentWager = Wager * 2
            end
            push!(BalResults, Balance)
        end

println()
#println("Gain/Loss after $GamesPlayed games: ", Balance - CashInHand,'$')
return BalResults
end


text:
using BayesianSumProductNetworks
using Test, Random
using AxisArrays, StatsFuns

priors = [NormalInverseGamma(0.0, 1.0, 1.0, 1.0), Beta()]
likelihoods = [Normal(), Bernoulli()]
sstats = [NormalInvGammaSufficientStats(), BetaSufficientStats()]
prior_sums_alpha = 1.0
prior_prods = 1.0
prior_leaves = 0.1
J = 2
K = 10
K_sum = 4
K_prod = 2

N = 100
D = 2

Random.seed!(1)
x = hcat(randn(N,1), rand(0:1, N,1))

# SPN
spn = templateRegion(prior_sums_alpha, prior_prods_alpha, prior_leaves,
    priors, likelihoods, sstats, N, D,
    K_sum, K_prod, J, K, 1, 2, root = true)

# Initialization
randomscopes!(spn, x)

@test scope(spn) == [1,2]

# Single ancestral sampling step
ancestralsampling!(spn, x)

@test sum(exp.(spn.logweights)) ≈ 1

# Single gibbs step
gibbssamplescopes!(spn, x)

@test scope(spn) == [1,2]

for child in children(spn)
    @test scope(child) == [1,2]
    scopes = map(n -> scope(n), children(child))
    @test length(reduce(vcat, scopes)) == 2
end

iterations = 1000

clusters = sort.([ [[1,2],[]], [[1],[2]] ])
lscopes = Vector{Vector}(undef, iterations)

for i in 1:iterations
    gibbssamplescopes!(spn, x)
    ancestralsampling!(spn, x)

    c = children(spn)[1]
    lscopes[i] = sort(map(n -> scope(n), children(c)))
end

cs = map(s -> findfirst(map(c -> c == s, clusters)), lscopes)

@test (sum(cs .== 1) / iterations) ≈ 0.5 atol=0.075
@test (sum(cs .== 2) / iterations) ≈ 0.5 atol=0.075

priors = [[NormalInverseGamma(0.0, 1.0, 1.0, 1.0), Gamma()], [Beta(), Gamma()]]
likelihoods = [Distribution[Normal(), Exponential()], Distribution[Bernoulli(), Poisson()]]
sstats = [[NormalInvGammaSufficientStats(), GammaSufficienStats()], [BetaSufficientStats(), GammaSufficienStats()]]

x[:,1] .-= minimum(x[:,1])

spn = templateRegion(prior_sums_alpha, prior_prods_alpha, prior_leaves,
    priors, likelihoods, sstats, N, D,
    K_sum, K_prod, J, K, 1, 2, root = true)

# Initialization
randomscopes!(spn, x)

@test scope(spn) == [1,2]

# Single ancestral sampling step
ancestralsampling!(spn, x)

@test sum(exp.(spn.logweights)) ≈ 1

# Single gibbs step
gibbssamplescopes!(spn, x)

@test scope(spn) == [1,2]

for child in children(spn)
    @test scope(child) == [1,2]
    scopes = map(n -> scope(n), children(child))
    @test length(reduce(vcat, scopes)) == 2
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
module ConstrainedFESpaces

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery

export ConstrainedFESpace

import Gridap: num_free_dofs
import Gridap: num_diri_dofs
import Gridap: diri_tags
import Gridap: apply_constraints
import Gridap: apply_constraints_rows
import Gridap: apply_constraints_cols
import Gridap: celldofids
import Gridap: interpolate_values
import Gridap: interpolate_diri_values
import Gridap: CellField
import Gridap: CellFieldForEval
import Gridap: CellBasis
import Gridap: Triangulation

using Base: @propagate_inbounds
import Base: size
import Base: IndexStyle
import Base: getindex

"""
A FESpace that it is constructed from another one plus some constraints
on the free values.
For the moment, we can only fix some free dofs. In the future, we will add
also linear constraints between free dofs.
"""
struct ConstrainedFESpace{D,Z,T,E} <: FESpace{D,Z,T}
  fespace::FESpace{D,Z,T}
  fixeddofs::Vector{Int}
  fixedvals::Vector{E}
  nfree::Int
  celldofs::IndexCellVector{Int}
  is_fixed::Vector{Bool}
  is_free::Vector{Bool}
  dof_to_new_dof::Vector{Int}
end

function ConstrainedFESpace(
  fespace::FESpace{D,Z,T},
  fixeddofs::Vector{Int},
  fixedvals::Vector{E}=zeros(dof_type(T),length(fixeddofs))) where {D,Z,T,E}

  @assert E == dof_type(T)

  celldofs, nfree, is_fixed, is_free, dof_to_new_dof =
    _setup_celldofs(fespace,fixeddofs)
  ConstrainedFESpace{D,Z,T,E}(
    fespace,fixeddofs,fixedvals,nfree,celldofs,
    is_fixed,is_free,dof_to_new_dof)
end

num_free_dofs(this::ConstrainedFESpace) = this.nfree

num_diri_dofs(this::ConstrainedFESpace) = num_diri_dofs(this.fespace)

diri_tags(this::ConstrainedFESpace) = diri_tags(this.fespace)

apply_constraints(
  this::ConstrainedFESpace,
  cellvec::CellVector,
  cellids::CellNumber) = apply_constraints(this.fespace,cellvec,cellids)

apply_constraints_rows(
  this::ConstrainedFESpace,
  cellmat::CellMatrix,
  cellids::CellNumber) = apply_constraints_rows(this.fespace,cellmat,cellids)

apply_constraints_cols(
  this::ConstrainedFESpace,
  cellmat::CellMatrix,
  cellids::CellNumber) = apply_constraints_cols(this.fespace,cellmat,cellids)

function celldofids(this::ConstrainedFESpace)
  this.celldofs
end

function interpolate_values(this::ConstrainedFESpace,fun::Function)
  freevals, dirivals = interpolate_values(this.fespace,fun)
  (freevals[this.is_free], dirivals)
end

function interpolate_values(this::ConstrainedFESpace{D,Z,T},val::T) where {D,Z,T}
  freevals, dirivals = interpolate_values(this.fespace,val)
  (freevals[this.is_free], dirivals)
end

function interpolate_diri_values(this::ConstrainedFESpace, funs::Vector{<:Function})
  interpolate_diri_values(this.fespace,funs)
end

function interpolate_diri_values(this::ConstrainedFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
  interpolate_diri_values(this.fespace,vals)
end

function CellField(
  this::ConstrainedFESpace{D,Z,T},
  free_dofs::AbstractVector{E},
  diri_dofs::AbstractVector{E}) where {D,Z,T,E}

  ndiri = length(diri_dofs)
  v = VectorOfTwoParts(this.dof_to_new_dof,free_dofs,this.fixedvals,ndiri)

  CellField(this.fespace,v,diri_dofs)
end

function CellFieldForEval(
  this::ConstrainedFESpace{D,Z,T},
  free_dofs::AbstractVector{E},
  diri_dofs::AbstractVector{E}) where {D,Z,T,E}

  ndiri = length(diri_dofs)
  v = VectorOfTwoParts(this.dof_to_new_dof,free_dofs,this.fixedvals,ndiri)

  CellFieldForEval(this.fespace,v,diri_dofs)
end

function CellBasis(this::ConstrainedFESpace{D,Z,T}) where {D,Z,T}
  CellBasis(this.fespace)
end

function Triangulation(this::ConstrainedFESpace{D,Z}) where {Z,D}
  Triangulation(this.fespace)
end

# Helpers

struct VectorOfTwoParts{T,VT1<:AbstractVector{T},VT2<:AbstractVector{T},I} <: AbstractVector{T}
  dof_to_xdof::Vector{I}
  pdof_to_val::VT1
  ndof_to_val::VT2
  offset::I
end

size(v::VectorOfTwoParts) = (length(v.dof_to_xdof),)

IndexStyle(::Type{VectorOfTwoParts{T,I}}) where {T,I} = IndexLinear()

@propagate_inbounds function getindex(v::VectorOfTwoParts,i::Integer)
  @inbounds xdof = v.dof_to_xdof[i]
  if xdof > 0
    @inbounds return v.pdof_to_val[xdof]
  elseif xdof< 0
    @inbounds return v.ndof_to_val[-xdof-v.offset]
  end
end

function _setup_celldofs(fespace,fixeddofs)

  celldofs = celldofids(fespace)
  nfree = num_free_dofs(fespace)
  ndiri = num_diri_dofs(fespace)

  nfixed = length(fixeddofs)

  nfree_new = nfree - nfixed

  dof_to_new_dof = zeros(Int,nfree)
  is_fixed = fill(false,nfree)
  is_fixed[fixeddofs] .= true
  is_free = is_fixed.==false

  dof_to_new_dof[is_free] .= 1:nfree_new
  dof_to_new_dof[is_fixed] .= (-(ndiri+1)):-1:(-(ndiri+nfixed))

  celldofs_new = CellVectorFromLocalToGlobalPosAndNeg(
    celldofs,dof_to_new_dof,collect(-1:-1:-ndiri))

  (celldofs_new, nfree_new, is_fixed, is_free, dof_to_new_dof)
end

end # module


text:
module CMPSKit

export virtualdim
export leftenv, rightenv, leftenv!, rightenv!, environments!
export leftgauge, rightgauge, leftgauge!, rightgauge!,
        leftcanonical, rightcanonical, leftcanonical!, rightcanonical!
export InfiniteCMPS, FiniteCMPS, LeftTransfer, RightTransfer
export Constant, FourierSeries, TaylorSeries, PiecewiseLinear
export fit, differentiate, integrate, localdot, domain, period, nummodes, density
export leftreducedoperator, rightreducedoperator, expval, gradient, groundstate
export groundstate2

export norm, dot, isapprox

export ∂, ∫, ψ̂, ∂ψ̂, ψ̂₁, ψ̂₂, ψ̂₃, ∂ψ̂₁, ∂ψ̂₂, ∂ψ̂₃

using LinearAlgebra
using KrylovKit
using KrylovKit: ConvergenceInfo
using OptimKit

scalartype(x::Any) = scalartype(typeof(x))
scalartype(T::Type{<:Number}) = T
scalartype(::Type{<:AbstractArray{T}}) where T = T

defaulttol(x::Any) = eps(real(float(one(scalartype(x)))))^(2/3)

isscalar(x) = false
isscalar(x::Number) = true

include("common/regularisedinverse.jl")
include("common/adjoint.jl")

include("functionspaces/generic.jl")
include("functionspaces/constant.jl")
include("functionspaces/fourierseries.jl")
include("functionspaces/taylorseries.jl")
include("functionspaces/piecewise.jl")
include("functionspaces/piecewiselinear.jl")

const MatrixFunction{T} = FunctionSpace{<:AbstractMatrix{T}}
const ScalarFunction = Union{<:Number, FunctionSpace{<:Number}}

include("operators/fieldoperators.jl")
include("operators/sumofterms.jl")
include("operators/hamiltonians.jl")

include("abstractcmps.jl")
include("transfer.jl")

include("finitecmps/finitecmps.jl")
include("finitecmps/transfer.jl")
include("finitecmps/environments.jl")
include("finitecmps/tangent.jl")
include("finitecmps/gradients.jl")
include("finitecmps/groundstate.jl")

include("infinitecmps/infinitecmps.jl")
include("infinitecmps/environments.jl")
include("infinitecmps/gauging.jl")
include("infinitecmps/gradients.jl")
include("infinitecmps/groundstate.jl")

# Deprecations
Base.@deprecate_binding ψ ψ̂
Base.@deprecate_binding ∂ψ ∂ψ̂

end


text:
module ProkarNC #prokaryotic auto-regulation
export iter, PAR, generate, runexamples, load
#using Winston

include("SDE.jl")
using .SDE
include("misc.jl")
import .SDE.sigma, .SDE.b, .SDE.a

######## data?
status = :nodata

####### VOU

immutable PAR <: MvPro
    th::Vector{Float64}
    d
    d2
    k::Float64   # number of copies of gene
    PAR(th) = new(th, 4, 8, 10.)
end
nthI(_::PAR) = 0
nthII(_::PAR) = 8
nth(_::PAR) = 8

const S =
   [ 0. 0  1  0  0  0 -1  0
     0  0  0  1 -2  2  0 -1
    -1  1  0  0  1 -1  0  0
    -1  1  0  0  0  0  0  0]
const St = S'

h(x, P::PAR) = max([max(x[3],0.)*x[4], P.k - x[4], x[4], x[1], 0.5*abs(x[2])*(x[2]-1),  x[3],  x[1],  x[2]], 0.)
thh(x, P::PAR) = P.th .* h(x,  P)
b(t,x, P::PAR) = S*thh(x, P)
sigma(t,x, P::PAR) = sqrt(max(thh(x, P),0.))' .* S
function a(t,x, P::PAR)
    s = sigma(t,x,P)
    s*s'
end

phi(t, x, P::PAR) = S .* h(x, P)'

#const U2 = 8.
#const I2 = -3.
#const U3 = 4.
#const U4 = 6.
#const I34 =0.

# 29.7586 5.02829 6.044613
# -52.47444 14.83634

#const U2 = 15.
#const I2 = -52.
#const U4 = 5.
#const U3 = 6.
#const I34 = 30.

const U2 = 18.
const I2 = -82.
const U4 = 5.
const U3 = 6.
const I34 = -32.

Bof(P) =
[   (-P.th[7])    0                   0                   P.th[3]
    P.th[4]     (-P.th[8]-U2*P.th[5]) 2P.th[6]            0
    0           0.5*U2*P.th[5]      (-P.th[6]-U4*P.th[1]) (-P.th[2] - U3*P.th[1])
    0           0                   -U4*P.th[1]           (-P.th[2] - U3*P.th[1])]

betaof(P) = [0, -I2*P.th[5],0.5I2*P.th[5] + P.k*P.th[2] - I34*P.th[1], P.k*P.th[2] - I34*P.th[1] ]

muof(P) = zeros(4)

####### Prior
import Base.clamp;
clamp(x::Vector,l::Vector, u::Vector) = map(i -> clamp(x[i], l[i], u[i]), 1:length(x))

function generate(Model_; thtrueI = [],
                          thtrueII = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1],
u = [7.5, 10., 5., 5.], Nobs = 51, N = 400000 + 1, Delta=1., SV=false)
    global Model = Model_
    th = [thtrueI, thtrueII]
    global Ptrue = Model(th)
    global model = string(Model)

    global M = Nobs - 1 # number of bridges
    global Ttotal = Delta*M      #total time span

    tag = "$model$(M)T$Ttotal"
    println(tag)
    d = Ptrue.d
    d2 = Ptrue.d2
    ############

    ############
    print("Generate X")
    ############

    srand(3)
    tt = linspace(0., Ttotal, N)
    W = sample(tt, Wiener(d2))
    xx = euler(u, W, Ptrue).yy
    println(".")

    global ttf = tt[1:(N-1)//M/100:end]
    global xxf = xx[:, 1:(N-1)//M/100:end]
    global ttd = tt[1:(N-1)//M:end]
    global xxd = xx[:, 1:(N-1)//M:end]

    global status = :generated

    if SV
        writecsv("xtrue$(tag).csv",[linspace(0, Ttotal, size(xx,2)) xx'])
        writecsv("xobs$(tag).csv",[linspace(0, Ttotal, size(xxd,2)) xxd'])
    end
end

function load(file="autoreg50fo", Model_=PAR; thtrueI = [],
            thtrueII = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1],
            Delta=1.)
    th = [thtrueI, thtrueII]

    const global Model = Model_
    const global Ptrue = Model(th)
    const global model = string(Model)

    const global xxd = readcsv(string(file,".csv"))'
    const global Nobs = size(xxd,2)
    const global M = Nobs - 1 # number of bridges
    const global Ttotal = Delta*M      #total time span
    const global ttd = linspace(0., Ttotal, Nobs)

    tag = "$model$file"
    println(tag)

    global status = :empir

end

function randq(th, delta, P)
        th .* exp(delta*(randn(length(th))))
end

function iter(K, n;             # numer of iterations, number of points for imputed bridges, including endpoints
proptype = :lin,                # type of proposals, :lin or :brown
TC = false,                     # apply timechange
PRECOMP = false,                # use tabulated Phim
SV = true, DRAW = false, VERBOSE = 1, REJNEG = true,
th1 = 0.05ones(8),
delta=0.05)
    NEEDX = true
    if status == :nodata error("Use generate(Model) or load(...) to generate or load observation process.") end

    tag = "NC$model$K$proptype$(M)x$(n)T$Ttotal"*(TC ? "TC":"")

    Nth = nth(Ptrue)
    NthII = nthII(Ptrue)
    NthI = nthI(Ptrue)
    d = Ptrue.d
    d2 = Ptrue.d2

## PRIOR
    pth(th) = 1.
    qth(th, thold) = 1.

    xi = 0.1*ones(NthI)

    ############
    VERBOSE > 0 && println("Generate bridges (using proposals $proptype): $tag.")
    ############
    srand(3)
    thetas = zeros(Nth, K)
    bb = zeros(Bool, M, K)
    bbp = zeros(M, K)
    aacc = zeros(Bool, K)
    aaccp = zeros(K)
    Td = zeros(M+1)

    llold = zeros(M)
    ll = zeros(M)
    ll2 = 0.

    laccth = 0.
    llmax = zeros(M)
    ll2max = 0
    rejected = 0
    rejectedth = 0

    TT = cell(M)
    TTU = cell(M)
    W = Array(CTPath{2},M)
    W° = Array(CTPath{2},M)
    Y = Array(CTPath{2},M)
    Y° = Array(CTPath{2},M)
    Z = Array(CTPath{2},M)
    Z° = Array(CTPath{2},M)
    U = Array(CTPath{2},M)
    U° = Array(CTPath{2},M)
    Phim = cell(n-1)
    Phim° = cell(n-1)

    # allocate arrays for bridges
    dt = Ttotal/M/n
    Tall = Ttotal/M
    TTUall = linspace(0.0, Tall, n)
    for m = 1:M
        Tmin = ttd[m]
        Tmax = ttd[m+1]
        TT[m] = linspace(Tmin, Tmax, n)
        TTU[m] = linspace(0.0, Tmax-Tmin, n)
        W[m] = MvPath(TT[m],d2)
        W°[m] = MvPath(TT[m],d2)
        Y[m] = MvPath(TT[m],d)
        Y°[m] = MvPath(TT[m],d)
        U[m] = MvPath(TTU[m],d)
        U°[m] = MvPath(TTU[m],d)
        Z[m] = MvPath(TTU[m],d2)
        Z°[m] = MvPath(TTU[m],d2)

    end

    P1 = Model(th1)
    th = copy(th1)
    th° = copy(th)
    println("th1: ", repr(th))
    qq = 1.0;

    srand(3)
    bbsumall = 0
    bbsum = zeros(M)
    aaccsum = 0
    mth = 0*th

    # conjugate posterior distribution N(mu, Si)
    mu = zeros(NthI)
    Si = zeros(NthI, NthI)

    for k = 1:K
        mth[:] = mth .+ th
        VERBOSE > 0 && k > 1 && print("k $k ", repr(round(mth./k,3)), " ",repr(round(th,3)))

        #### generate bridges
        P = Model(th)
        if PRECOMP
            for i = 1:n-1
                Phim[i] = Phims(TTUall[i], Tall, Bof(P))
            end
        end
        for m = 1:M

            u = xxd[:, m]
            v = xxd[:, m+1]

            Tmin = ttd[m]
            Tmax = ttd[m+1]
            T = Tmax-Tmin

            if proptype == :brown || proptype == :brownbrown
                Pt = MvAffPro(muof(P), sigma(Tmax,v, P))
                #Pt  = MvLinPro(-1E-7, 0.0, sigma(Tmax,v, P))
            elseif proptype == :lin || proptype == :linlin
                Pt = MvLinPro(Bof(P),betaof(P), sigma(Tmax,v,P))
            end

            posi = false
            while !posi

                if TC && PRECOMP
                    resample!(Z°[m], Wiener(d2))
                    u0 = uofx(0., u, T, v, Pt, Phim[1])
                    eulerU!(U°[m], u0, Z°[m], Tmin, T, v, Pt, P, Phim)
    NEEDX &&        XofU!(Y°[m], U°[m], Tmin, T, v, Pt)
                    ll[m] = llikeliU(U°[m], Tmin, T, v, Pt, P, Phim)
                    llold[m] = llikeliU(U[m], Tmin, T, v, Pt, P, Phim)
                elseif TC
                    resample!(Z°[m], Wiener(d2))
                    u0 = uofx(0., u, T, v, Pt)
                    eulerU!(U°[m], u0, Z°[m], Tmin, T, v, Pt, P)
    NEEDX &&        XofU!(Y°[m], U°[m], Tmin, T, v, Pt)
                    ll[m] = llikeliU(U°[m], Tmin, T, v, Pt, P)
					llold[m] = llikeliU(U[m], Tmin, T, v, Pt, P)
                else
                    resample!(W°[m], Wiener(d2))
                    guidedeuler!(Y°[m], u, W°[m], Tmax, v, Pt, P)
                    ll[m] = llikeliXcirc(Y°[m], Pt, P)
                    llold[m] = llikeliXcirc(Y[m], Pt, P)
                end
                posi = true
                if REJNEG && any(Y°[m].yy .<= 0)
                    posi = false
                    print("\$")
                    rejected += 1
                end

            end

            llmax = map(max, ll, llmax)
            if(k == 1)
                    Z[m].yy[:] = Z°[m].yy
                    U[m].yy[:] = U°[m].yy
                    Y[m].yy[:] = Y°[m].yy
                    Y[m].tt[:] = Y°[m].tt
                    W[m].yy[:] = W°[m].yy
                    llold[m] = ll[m]
            end

            acc = min(1.0, exp(ll[m]-llold[m]))
            VERBOSE > 1 && println("\t acc ", round(acc, 3), " ", round(llold[m],4), " ", round(ll[m],4))
            bbp[m, k] = acc
            bbsum[m] += acc
            bbsumall += acc
            bb[m, k] = false
            if rand() <= acc
                bb[m, k] = true
	            Y[m], Y°[m] = Y°[m], Y[m]
                W[m], W°[m] = W°[m], W[m]
	            Z[m], Z°[m] = Z°[m], Z[m]
                U[m], U°[m] = U°[m], U[m]
            end
        end #for m

        if DRAW
             I1 = 1 + k % 40
             I2 = I1 + 9
             p = FramedPlot()
             setattr(p, "xrange", (I1, I2))
             setattr(p, "yrange", (-3,20))

             for j in 1:d
                 add(p, Points(ttd[I1:I2], xxd[j, I1:I2+1][:], "color","red"))
                 for i in I1:I2
                     add(p, Curve(Y[i].tt[:],Y[i].yy[j, :][:] , "color","black", "linewidth", 0.5))
                 end
             end
             display(p)
        end
        if false
             C1 = 3
             C2 = 4
             if k % 4 == 0; C1 = 1; C2 = 4; end
             I1 = 1 + k % 40
             I2 = I1 + 9
             p = FramedPlot()
             const Rs = [(-1, 15),(-1, 30),(-5, 15),(-5, 15)]
             setattr(p, "xrange", Rs[C1])
             setattr(p, "yrange", Rs[C2])
             add(p, Points(xxd[C1, I1:I2+1][:], xxd[C2, I1:I2+1][:], "color","red"))
             for i in I1:I2
                 add(p, Curve(Y[i].yy[C1, :][:],Y[i].yy[C2, :][:] , "color","black", "linewidth", 0.5))
             end
             display(p)
        end

		VERBOSE > 0 && k > 1 && print(" ", round(100*bbsumall/k/M,1),"(min ", round(100*minimum(bbsum)/k,1))
		VERBOSE > 0 && k > 1 && @printf(" ll %.1f)",round(maximum(ll),1))

		P = Model(th)

        posi = false
        while !posi
            posi = true
            th° = randq(th, delta, P)
            th° = min(1100.,max(0.00091, th°))

			P° = Model(th°)
		    VERBOSE > 0 && print(" th° ", repr(round(th°,3)))

            if PRECOMP
                for i = 1:n-1
                    Phim°[i] = Phims(TTUall[i], Tall, Bof(P°))
                end
            end

            #### update sigma and theta
            laccth = 0.

            for m in 1:M
                u = xxd[:, m]
                v = xxd[:, m+1]

                Tmin = ttd[m]
                Tmax = ttd[m+1]
                T = Tmax-Tmin

                if proptype == :brown || proptype == :brownbrown
                    Pt = MvAffPro(muof(P), sigma(Tmax,v, P))
                    Pt° = MvAffPro(muof(P°),  sigma(Tmax,v, P°))
                elseif proptype == :lin || proptype == :linlin
                    Pt = MvLinPro(Bof(P),betaof(P), sigma(Tmax,v,P))
                    Pt° = MvLinPro(Bof(P°),betaof(P°), sigma(Tmax,v, P°))
                end

                # compute acceptance probability of change in theta with equivalent Innovations

                ### bridge Y° (U) using original W (Z)
                if TC && PRECOMP
                    u0 = uofx(0., u, T, v, Pt°, Phim°[1])
                    eulerU!(U°[m], u0, Z[m], Tmin, T, v, Pt°, P°, Phim°)
    NEEDX &&        XofU!(Y°[m], U°[m], Tmin, T, v, Pt°)
                    m1 = +lp(Tmin, u, Tmax, v, Pt°)
                    m2 = -lp(Tmin, u, Tmax, v, Pt)
                    m3 = +llikeliU(U°[m], Tmin, T, v, Pt°, P°, Phim°)
                    m4 = -llikeliU(U[m], Tmin, T, v, Pt, P, Phim)
                    laccth += m1 + m2 + m3 + m4
                elseif TC
                    u0 = uofx(0., u, T, v, Pt°)
                    eulerU!(U°[m], u0, Z[m], Tmin, T, v, Pt°, P°)
    NEEDX &&        XofU!(Y°[m], U°[m], Tmin, T, v, Pt°)
                    m1 = +lp(Tmin, u, Tmax, v, Pt°)
                    m2 = -lp(Tmin, u, Tmax, v, Pt)
                    m3 = +llikeliU(U°[m], Tmin, T, v, Pt°, P°)
                    m4 = -llikeliU(U[m], Tmin, T, v, Pt, P)
                    laccth += m1 + m2 + m3 + m4
                else
                    guidedeuler!(Y°[m], u, W[m], Tmax, v, Pt°, P°)
                    m1 = +lp(Tmin, u, Tmax, v, Pt°)
                    m2 = -lp(Tmin, u, Tmax, v, Pt)
                    m3 = +llikeliXcirc(Y°[m], Pt°, P°)
                    m4 = -llikeliXcirc(Y[m], Pt, P)
                    laccth += m1 + m2 + m3 + m4
                end
                if REJNEG && any(Y°[m].yy .<= 0)
                        posi = false
                        print("\$")
                        rejectedth += 1
                end

            end #for m

        end #while !posi

        ll2max = max(laccth, ll2max)
        accth = min(1., exp(laccth) * pth(th°)/pth(th)*qth(th, th°)/qth(th°, th))
        aacc[k] = false
        aaccp[k] = accth
        VERBOSE > 0 && print(" Pacc $(repr(round(accth,3))) avg ", round(100aaccsum/k,2))

        if rand() <= accth
            aacc[k] = true
            aaccsum += 1

            th[:] = th°
            for m = 1:M
                # W and Z do not change
                Y[m], Y°[m] = Y°[m], Y[m]
                U[m], U°[m] = U°[m], U[m]
            end

        end
		VERBOSE > 0 && @printf(" ll %.1f %s", round(laccth,1), aacc[k] ? "acc" : "")

		thetas[:, k] = th
		println()

    end # for K

    if SV
        if !isdir(tag)
             mkdir(tag)
        end
        writecsv("$tag/thetas$tag.csv", thetas')
        accs = [aacc aaccp bb' bbp']
        writecsv("$tag/accs$tag.csv", accs)
        open("$tag/acc$tag.txt", "w") do f
            println(f,"estimate: $(repr(mean(thetas,2)[:]))")
            println(f,"bridge acc ",  repr(round(100*mean(1.0bb),2)), "%, avg prob  ", repr(round(100*mean(bbp),2)) )
            println(f, "llmax ", maximum(llmax), "rejected neg. path", rejected)
            println(f,"sigma acc ",  round(100*sum(aacc)/K,2), "% avg prob ", round(100*sum(aaccp)/K, 2), "% " )
        end
    end

    println("estimate: $(repr(mean(thetas,2)[:]))")
    println("bridge acc ",  round(100*mean(1.0bb,2),2), "%, avg prob  ", repr(round(100*mean(bbp, 2),2)) )
    println("sigma acc ",  round(100*sum(aacc)/K,2), "% avg prob ", round(100*sum(aaccp)/K, 2), "% " )
    return thetas, accs
end

function runexamples(k=100000)
    load("autoreg50fo", PAR)
    for n in [20,50,10], t in [:lin], tc in [true]
        println("$n $t $tc")
        iter(k, n; proptype=t, DRAW=false, TC=tc, PRECOMP=true, delta=.12, REJNEG=true)
    end
end

end #module


text:
function get_current_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getCurrentState", ODEStateAndDerivative, ())
end

function get_global_current_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getGlobalCurrentState", ODEStateAndDerivative, ())
end

function get_global_previous_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getGlobalPreviousState", ODEStateAndDerivative, ())
end

function get_interpolated_state(obj::AbstractODEStateInterpolator, arg0::jdouble)
    return jcall(obj, "getInterpolatedState", ODEStateAndDerivative, (jdouble,), arg0)
end

function get_previous_state(obj::AbstractODEStateInterpolator)
    return jcall(obj, "getPreviousState", ODEStateAndDerivative, ())
end

function is_current_state_interpolated(obj::AbstractODEStateInterpolator)
    return jcall(obj, "isCurrentStateInterpolated", jboolean, ())
end

function is_forward(obj::AbstractODEStateInterpolator)
    return jcall(obj, "isForward", jboolean, ())
end

function is_previous_state_interpolated(obj::AbstractODEStateInterpolator)
    return jcall(obj, "isPreviousStateInterpolated", jboolean, ())
end

function restrict_step(obj::AbstractODEStateInterpolator, arg0::ODEStateAndDerivative, arg1::ODEStateAndDerivative)
    return jcall(obj, "restrictStep", AbstractODEStateInterpolator, (ODEStateAndDerivative, ODEStateAndDerivative), arg0, arg1)
end


text:
import Adapt: adapt, adapt_storage
using  LinearAlgebra: Cholesky
using Zygote: IdSet
import Functors: Functors, @functor, functor, fmap, isleaf
using SparseArrays: AbstractSparseArray

trainable(m) = functor(m)[1]

"""
    testmode!(m, mode = true)

Set a layer or model's test mode (see below).
Using `:auto` mode will treat any gradient computation as training.

_Note_: if you manually set a model into test mode, you need to manually place
it back into train mode during training phase.

Possible values include:
- `false` for training
- `true` for testing
- `:auto` or `nothing` for Flux to detect the mode automatically
"""
testmode!(m, mode = true) = (foreach(x -> testmode!(x, mode), trainable(m)); m)

"""
    trainmode!(m, mode = true)

Set a layer of model's train mode (see below).
Symmetric to [`testmode!`](@ref) (i.e. `trainmode!(m, mode) == testmode!(m, !mode)`).

_Note_: if you manually set a model into train mode, you need to manually place
it into test mode during testing phase.

Possible values include:
- `true` for training
- `false` for testing
- `:auto` or `nothing` for Flux to detect the mode automatically
"""
trainmode!(m, mode = true) = mode isa Bool ? testmode!(m, !mode) : testmode!(m, mode)

params!(p::Params, x::AbstractArray{<:Number}, seen = IdSet()) = push!(p, x)

function params!(p::Params, x, seen = IdSet())
  x in seen && return
  push!(seen, x)
  for child in trainable(x)
    params!(p, child, seen)
  end
end

"""
    params(model)
    params(layers...)

Given a model or specific layers from a model, create a `Params` object pointing to its trainable parameters.

This can be used with the `gradient` function, see [Taking Gradients](@ref), or as input to the [`Flux.train!`](@ref Flux.train!) function.

The behaviour of `params` on custom types can be customized using [`Functor.@functor`](@ref) or [`Flux.trainable`](@ref).

# Examples
```jldoctest
julia> using Flux: params

julia> params(Chain(Dense(ones(2,3)), softmax))  # unpacks Flux models
Params([[1.0 1.0 1.0; 1.0 1.0 1.0], [0.0, 0.0]])

julia> bn = BatchNorm(2, relu)
BatchNorm(2, relu)  # 4 parameters, plus 4 non-trainable

julia> params(bn)  # only the trainable parameters
Params([Float32[0.0, 0.0], Float32[1.0, 1.0]])

julia> params([1, 2, 3], [4])  # one or more arrays of numbers
Params([[1, 2, 3], [4]])

julia> params([[1, 2, 3], [4]])  # unpacks array of arrays
Params([[1, 2, 3], [4]])

julia> params(1, [2 2], (alpha=[3,3,3], beta=Ref(4), gamma=sin))  # ignores scalars, unpacks NamedTuples
Params([[2 2], [3, 3, 3]])
```
"""
function params(m...)
  ps = Params()
  params!(ps, m)
  return ps
end

function loadparams!(m, xs)
  for (p, x) in zip(params(m), xs)
    size(p) == size(x) ||
      error("Expected param size $(size(p)), got $(size(x))")
    copyto!(p, x)
  end
end

struct FluxCUDAAdaptor end
adapt_storage(to::FluxCUDAAdaptor, x) = CUDA.cu(x)
adapt_storage(to::FluxCUDAAdaptor, x::Zygote.FillArrays.AbstractFill) = CUDA.cu(collect(x))
if VERSION >= v"1.7"
  adapt_storage(to::FluxCUDAAdaptor, x::Random.TaskLocalRNG) = CUDA.default_rng()
else
  adapt_storage(to::FluxCUDAAdaptor, x::Random._GLOBAL_RNG) = CUDA.default_rng()
end
adapt_storage(to::FluxCUDAAdaptor, x::CUDA.RNG) = x
adapt_storage(to::FluxCUDAAdaptor, x::AbstractRNG) =
  error("Cannot map RNG of type $(typeof(x)) to GPU. GPU execution only supports Random.default_rng().")

# TODO: figure out the correct design for OneElement
adapt_storage(to::FluxCUDAAdaptor, x::Zygote.OneElement) = CUDA.cu(collect(x))

struct FluxCPUAdaptor end

# define rules for handling structured arrays
adapt_storage(to::FluxCPUAdaptor, x::AbstractArray) = adapt(Array, x)
adapt_storage(to::FluxCPUAdaptor, x::AbstractRange) = x
adapt_storage(to::FluxCPUAdaptor, x::Zygote.FillArrays.AbstractFill) = x
adapt_storage(to::FluxCPUAdaptor, x::T) where T <: CUDA.CUSPARSE.CUDA.CUSPARSE.AbstractCuSparseMatrix = adapt(Array, x)
adapt_storage(to::FluxCPUAdaptor, x::Zygote.OneElement) = x
adapt_storage(to::FluxCPUAdaptor, x::AbstractSparseArray) = x
adapt_storage(to::FluxCPUAdaptor, x::CUDA.RNG) = Random.default_rng()
adapt_storage(to::FluxCPUAdaptor, x::AbstractRNG) = x

function ChainRulesCore.rrule(::typeof(Array), x::CUDA.CuArray)
  Array(x), d -> (NoTangent(), CUDA.cu(d),)
end

function ChainRulesCore.rrule(::typeof(Adapt.adapt_storage), to::FluxCPUAdaptor, x::CUDA.AbstractGPUArray)
  adapt_storage(to, x), d -> (NoTangent(), NoTangent(), adapt_storage(FluxCUDAAdaptor(), d),)
end

# CPU/GPU movement conveniences

"""
    cpu(m)

Moves `m` onto the CPU, the opposite of [`gpu`](@ref).
Recurses into structs marked [`@functor`](@ref).

```julia-repl
julia> m = Dense(1,2)
Dense(1, 2)

julia> m_gpu = gpu(m)
Dense(1, 2)

julia> typeof(m_gpu.W)
CuArray{Float32, 2}

julia> m_cpu = cpu(m_gpu)
Dense(1, 2)

julia> typeof(m_cpu.W)
Matrix{Float32}
```
"""
cpu(x) = fmap(x -> adapt(FluxCPUAdaptor(), x), x)

_isbitsarray(::AbstractArray{<:Number}) = true
_isbitsarray(::AbstractArray{T}) where T = isbitstype(T)
_isbitsarray(x) = false

_isleaf(::AbstractRNG) = true
_isleaf(x) = _isbitsarray(x) || Functors.isleaf(x)

"""
    gpu(x)

Moves `m` to the current GPU device, if available. It is a no-op otherwise.
See the [CUDA.jl docs](https://juliagpu.github.io/CUDA.jl/stable/usage/multigpu/)
to help identify the current device.

This works for functions, and any struct marked with [`@functor`](@ref).

```julia-repl
julia> m = Dense(1,2)
Dense(1, 2)

julia> typeof(m.W)
Matrix{Float32}

julia> m_gpu = gpu(m)
Dense(1, 2)

julia> typeof(m_gpu.W) # notice the type of the array changed to a CuArray
CuArray{Float32, 2}
```
"""
function gpu(x)
  check_use_cuda()
  use_cuda[] ? fmap(x -> Adapt.adapt(FluxCUDAAdaptor(), x), x; exclude = _isleaf) : x
end

function check_use_cuda()
  if use_cuda[] === nothing
    use_cuda[] = CUDA.functional()
    if use_cuda[] && !CUDA.has_cudnn()
      @warn "CUDA.jl found cuda, but did not find libcudnn. Some functionality will not be available."
    end
    if !(use_cuda[])
      @info """The GPU function is being called but the GPU is not accessible.
               Defaulting back to the CPU. (No action is required if you want to run on the CPU).""" maxlog=1
    end
  end
end
ChainRulesCore.@non_differentiable check_use_cuda()

# Precision

adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Real}) = convert.(T, xs) # piracy

paramtype(T::Type{<:Real}, m) = fmap(x -> adapt(T, x), m)

"""
    f32(m)

Convert the `eltype` of model's parameters to `Float32`.
"""
f32(m) = paramtype(Float32, m)

"""
    f64(m)

Convert the `eltype` of model's parameters to `Float64`.
"""
f64(m) = paramtype(Float64, m)

# Functors for certain Julia data structures
@functor Cholesky
trainable(c::Cholesky) = ()


text:
include("imgui_widgets.jl")

using BangBang

function render(s::PupperSim.mjSim, w::GLFW.Window)
    # Update scene
    mjv_updateScene(s.m, s.d, s.vopt, s.pert, s.cam, MJCore.mjCAT_ALL, s.scn)

    # acquire head mounted camera image
    s.showcam && get_cameraview_image(s)    # must occur before main screen render

    # Render
    width, height = GLFW.GetFramebufferSize(w)

    #println("shadowClip: $(s.con.x.shadowClip), shadowScale: $(s.con.x.shadowScale), shadowSize: $(s.con.x.shadowSize)")

    elapsed_time = @elapsed mjr_render(mjrRect(0,0,width,height), s.scn, s.con)
    #println(round(elapsed_time * 1000; digits=3), " ms")

    #println(s.scn[].flags[1])  # shadows on (1) or off (0)
    if s.scn[].flags[1] == 1
        if elapsed_time > 0.01  # Rendering time is greater than 10 ms
            if s.m.vis.quality.shadowsize > 64
                @set!! s.m.vis.quality.shadowsize = s.m.vis.quality.shadowsize / 2
                mjr_makeContext(s.m, s.con, FONTSCALE)
                #println("Decreased shadowsize: $(s.m.vis.quality.shadowsize)")
            else
                @info("Turning off shadows to speed up rendering")
                flags = MVector(s.scn[].flags)
                flags[1] = 0
                s.scn[].flags = flags
            end
        elseif elapsed_time < 0.002 && s.m.vis.quality.shadowsize < 16384
            @set!! s.m.vis.quality.shadowsize = s.m.vis.quality.shadowsize * 2
            mjr_makeContext(s.m, s.con, FONTSCALE)  # model specific setup
            #println("Increased shadowsize: $(s.m.vis.quality.shadowsize)")
        end
    end

    #mjui_render(s.ui1, s.uistate, s.con)

    # display head mounted camera image
    s.showcam && render_camera_view(s)  # must occur after main screen render

    # If recording, record frame this frame
    s.record === nothing || record_video_frame(s, width, height)

    # Swap front and back buffers
    GLFW.SwapBuffers(w)
end


text:
s = open("input.txt") do file
    read(file, String);
end

groups = split(s,"\n\n");

function GetQuestionCount(group)
    people = split(strip(group),"\n");

    questions = Set(people[1]);

    for i = 2:length(people)
        intersect!(questions, Set(people[i]));
    end

    return length(questions);
end

sum(map(GetQuestionCount, groups))


text:
include("build_Zlib.v1.2.11.jl")
include("build_XML2.v2.9.9.jl")