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

mutable struct ManagementPolicyAction <: SwaggerModel
    baseBlob::Any # spec type: Union{ Nothing, ManagementPolicyBaseBlob } # spec name: baseBlob
    snapshot::Any # spec type: Union{ Nothing, ManagementPolicySnapShot } # spec name: snapshot

    function ManagementPolicyAction(;baseBlob=nothing, snapshot=nothing)
        o = new()
        validate_property(ManagementPolicyAction, Symbol("baseBlob"), baseBlob)
        setfield!(o, Symbol("baseBlob"), baseBlob)
        validate_property(ManagementPolicyAction, Symbol("snapshot"), snapshot)
        setfield!(o, Symbol("snapshot"), snapshot)
        o
    end
end # type ManagementPolicyAction

const _property_map_ManagementPolicyAction = Dict{Symbol,Symbol}(Symbol("baseBlob")=>Symbol("baseBlob"), Symbol("snapshot")=>Symbol("snapshot"))
const _property_types_ManagementPolicyAction = Dict{Symbol,String}(Symbol("baseBlob")=>"ManagementPolicyBaseBlob", Symbol("snapshot")=>"ManagementPolicySnapShot")
Base.propertynames(::Type{ ManagementPolicyAction }) = collect(keys(_property_map_ManagementPolicyAction))
Swagger.property_type(::Type{ ManagementPolicyAction }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagementPolicyAction[name]))}
Swagger.field_name(::Type{ ManagementPolicyAction }, property_name::Symbol) =  _property_map_ManagementPolicyAction[property_name]

function check_required(o::ManagementPolicyAction)
    true
end

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


text:
module OpenColor
  using ColorTypes

  export OCGRAY, OCRED, OCPINK, OCGRAPE,
    OCVIOLET, OCINDIGO, OCBLUE, OCCYAN,
    OCTEAL, OCGREEN, OCYELLOW, OCORANGE

  # Use Dict instead of NamedTuple because of zero-based index
  const OCGRAY = Base.Dict(0:9 .=> [
    RGB(248/255, 249/255, 250/255),
    RGB(241/255, 243/255, 245/255),
    RGB(233/255, 236/255, 239/255),
    RGB(222/255, 226/255, 230/255),
    RGB(206/255, 212/255, 218/255),
    RGB(173/255, 181/255, 189/255),
    RGB(134/255, 142/255, 150/255),
    RGB(73/255, 80/255, 87/255),
    RGB(52/255, 58/255, 64/255),
    RGB(33/255, 37/255, 41/255)])

  const OCRED = Base.Dict(0:9 .=> [
    RGB(255/255, 245/255, 245/255),
    RGB(255/255, 227/255, 227/255),
    RGB(255/255, 201/255, 201/255),
    RGB(255/255, 168/255, 168/255),
    RGB(255/255, 135/255, 135/255),
    RGB(255/255, 107/255, 107/255),
    RGB(250/255, 82/255, 82/255),
    RGB(240/255, 62/255, 62/255),
    RGB(224/255, 49/255, 49/255),
    RGB(201/255, 42/255, 42/255)])

  const OCPINK = Base.Dict(0:9 .=> [
    RGB(255/255, 240/255, 246/255),
    RGB(255/255, 222/255, 235/255),
    RGB(252/255, 194/255, 215/255),
    RGB(250/255, 162/255, 193/255),
    RGB(247/255, 131/255, 172/255),
    RGB(240/255, 101/255, 149/255),
    RGB(230/255, 73/255, 128/255),
    RGB(214/255, 51/255, 108/255),
    RGB(194/255, 37/255, 92/255),
    RGB(166/255, 30/255, 77/255)])

  const OCGRAPE = Base.Dict(0:9 .=> [
    RGB(248/255, 240/255, 252/255),
    RGB(243/255, 217/255, 250/255),
    RGB(238/255, 190/255, 250/255),
    RGB(229/255, 153/255, 247/255),
    RGB(218/255, 119/255, 242/255),
    RGB(204/255, 93/255, 232/255),
    RGB(190/255, 75/255, 219/255),
    RGB(174/255, 62/255, 201/255),
    RGB(156/255, 54/255, 181/255),
    RGB(134/255, 46/255, 156/255)])

  const OCVIOLET = Base.Dict(0:9 .=> [
    RGB(243/255, 240/255, 255/255),
    RGB(229/255, 219/255, 255/255),
    RGB(208/255, 191/255, 255/255),
    RGB(177/255, 151/255, 252/255),
    RGB(151/255, 117/255, 250/255),
    RGB(132/255, 94/255, 247/255),
    RGB(121/255, 80/255, 242/255),
    RGB(112/255, 72/255, 232/255),
    RGB(103/255, 65/255, 217/255),
    RGB(95/255, 61/255, 196/255)])

  const OCINDIGO = Base.Dict(0:9 .=> [
    RGB(237/255, 242/255, 255/255),
    RGB(219/255, 228/255, 255/255),
    RGB(186/255, 200/255, 255/255),
    RGB(145/255, 167/255, 255/255),
    RGB(116/255, 143/255, 252/255),
    RGB(92/255, 124/255, 250/255),
    RGB(76/255, 110/255, 245/255),
    RGB(66/255, 99/255, 235/255),
    RGB(59/255, 91/255, 219/255),
    RGB(54/255, 79/255, 199/255)])

  const OCBLUE = Base.Dict(0:9 .=> [
    RGB(231/255, 245/255, 255/255),
    RGB(208/255, 235/255, 255/255),
    RGB(165/255, 216/255, 255/255),
    RGB(116/255, 192/255, 252/255),
    RGB(77/255, 171/255, 247/255),
    RGB(51/255, 154/255, 240/255),
    RGB(34/255, 139/255, 230/255),
    RGB(28/255, 126/255, 214/255),
    RGB(25/255, 113/255, 194/255),
    RGB(24/255, 100/255, 171/255)])

  const OCCYAN = Base.Dict(0:9 .=> [
    RGB(227/255, 250/255, 252/255),
    RGB(197/255, 246/255, 250/255),
    RGB(153/255, 233/255, 242/255),
    RGB(102/255, 217/255, 232/255),
    RGB(59/255, 201/255, 219/255),
    RGB(34/255, 184/255, 207/255),
    RGB(21/255, 170/255, 191/255),
    RGB(16/255, 152/255, 173/255),
    RGB(12/255, 133/255, 153/255),
    RGB(11/255, 114/255, 133/255)])

  const OCTEAL = Base.Dict(0:9 .=> [
    RGB(230/255, 252/255, 245/255),
    RGB(195/255, 250/255, 232/255),
    RGB(150/255, 242/255, 215/255),
    RGB(99/255, 230/255, 190/255),
    RGB(56/255, 217/255, 169/255),
    RGB(32/255, 201/255, 151/255),
    RGB(18/255, 184/255, 134/255),
    RGB(12/255, 166/255, 120/255),
    RGB(9/255, 146/255, 104/255),
    RGB(8/255, 127/255, 91/255)])

  const OCGREEN = Base.Dict(0:9 .=> [
    RGB(235/255, 251/255, 238/255),
    RGB(211/255, 249/255, 216/255),
    RGB(178/255, 242/255, 187/255),
    RGB(140/255, 233/255, 154/255),
    RGB(105/255, 219/255, 124/255),
    RGB(81/255, 207/255, 102/255),
    RGB(64/255, 192/255, 87/255),
    RGB(55/255, 178/255, 77/255),
    RGB(47/255, 158/255, 68/255),
    RGB(43/255, 138/255, 62/255)])

  const OCLIME = Base.Dict(0:9 .=> [
    RGB(244/255, 252/255, 227/255),
    RGB(233/255, 250/255, 200/255),
    RGB(216/255, 245/255, 162/255),
    RGB(192/255, 235/255, 117/255),
    RGB(169/255, 227/255, 75/255),
    RGB(148/255, 216/255, 45/255),
    RGB(130/255, 201/255, 30/255),
    RGB(116/255, 184/255, 22/255),
    RGB(102/255, 168/255, 15/255),
    RGB(92/255, 148/255, 13/255)])

  const OCYELLOW = Base.Dict(0:9 .=> [
    RGB(255/255, 249/255, 219/255),
    RGB(255/255, 243/255, 191/255),
    RGB(255/255, 236/255, 153/255),
    RGB(255/255, 224/255, 102/255),
    RGB(255/255, 212/255, 59/255),
    RGB(252/255, 196/255, 25/255),
    RGB(250/255, 176/255, 5/255),
    RGB(245/255, 159/255, 0/255),
    RGB(240/255, 140/255, 0/255),
    RGB(230/255, 119/255, 0/255)])

  const OCORANGE = Base.Dict(0:9 .=> [
    RGB(255/255, 244/255, 230/255),
    RGB(255/255, 232/255, 204/255),
    RGB(255/255, 216/255, 168/255),
    RGB(255/255, 192/255, 120/255),
    RGB(255/255, 169/255, 77/255),
    RGB(255/255, 146/255, 43/255),
    RGB(253/255, 126/255, 20/255),
    RGB(247/255, 103/255, 7/255),
    RGB(232/255, 89/255, 12/255),
    RGB(217/255, 72/255, 15/255)])
end

ext: jl


text:
using Documenter, DocumenterMarkdown, CyclotomicNumbers

makedocs(sitename="CyclotomicNumbers documentation",format=Markdown(),modules=[CyclotomicNumbers])


text:
using SpecialFunctions: loggamma
using PyPlot

include("assess.jl")

α_vals = Float64[ n/50 for n = 10:2:50 ]
β = 1.0
tol = 1e-16
x = range(1, 50, length=50)
Nmin = Array{Float64}(undef, 50)

figure(1)

for α in α_vals
    for k = 1:50
        Nmin[k], term = assess(α, β, x[k], tol)
        if Nmin[k] < 0
            Nmin[k] = NaN
        end
    end
    plot(x, Nmin, label="α = $α")
end
title("tol = $tol")
#legend()
xlabel("x")
ylabel("N")
grid(true)


text:
@testset "Minimal Dominating Set" begin

    g0 = SimpleGraph(0)
    for g in testgraphs(g0)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        @test isempty(d)
    end

    g1 = SimpleGraph(1)
    for g in testgraphs(g1)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        @test (d == [1,])
    end

    add_edge!(g1, 1, 1)
    for g in testgraphs(g1)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        @test (d == [1,])
    end

    g3 = StarGraph(5)
    for g in testgraphs(g3)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        @test (length(d)== 1 || (length(d)== 4 && minimum(d) > 1 ))
    end

    g4 = CompleteGraph(5)
    for g in testgraphs(g4)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        @test length(d)== 1 #Exactly one vertex
    end

    g5 = PathGraph(4)
    for g in testgraphs(g5)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        sort!(d)
        @test (d == [1, 2, 4] || d == [1, 3] || d == [1, 4] || d == [2, 3] || d == [2, 4])
    end

    add_edge!(g5, 2, 2)
    add_edge!(g5, 3, 3)
    for g in testgraphs(g5)
        d = @inferred(dominating_set(g, MinimalDominatingSet(); seed=3))
        sort!(d)
        @test (d == [1, 2, 4] || d == [1, 3] || d == [1, 4] || d == [2, 3] || d == [2, 4])
    end
end


text:
#=
This is an example of how to use a custom user type of with ODE module.

An additional example can be found in `test/inverface-tests.jl`

The ODE considered in this example is the simple oscilator

  a' = -b
  b' =  a

with initial condtion a(t=0) = 0 and b(t=0) = 0.1.
=#

using ODE
using PyPlot
using Compat

const DataFloat = Float64
mutable struct Vec2
  a::DataFloat
  b::DataFloat
end

# This is the minimal set of function required on the type inorder to work with
# the ODE module
Base.:/(x::Vec2, y::Real) = Vec2(x.a/y, x.b/y)
Base.:*(y::Real, x::Vec2) = Vec2(y * x.a, y * x.b)
Base.:*(x::Vec2, y::Real) = y*x
Base.:.*(y::Real, x::Vec2) = y*x
Base.:+(x::Vec2, y::Vec2) = Vec2(x.a + y.a, x.b + y.b)
Base.:-(x::Vec2, y::Vec2) = Vec2(x.a - y.a, x.b - y.b)
Base.norm(x::Vec2) = sqrt(x.a^2 + x.b^2)
Base.zero(x::Type{Vec2}) = Vec2(zero(DataFloat), zero(DataFloat))
ODE.isoutofdomain(x::Vec2) = isnan(x.a) || isnan(x.b)

# RHS function
f(t,y) = Vec2(-y.b, y.a)

# Initial condtions
start = Vec2(0.0, 0.1)

# Time vector going from 0 to 2*PI in 0.01 increments
time = 0:0.1:4*pi

# Solve the ODE
t, y = ode45(f, start, time)

# Plot the solution
a = map(y -> y.a, y)
b = map(y -> y.b, y)
plot(t, a, label="a(t)")
plot(t, b, label="b(t)")
legend()


text:
using DiscreteChoice
using Optim
using ForwardDiff
using Random, StatsBase

#################
# Simulate data
#################

# Simulate J=100 case

Random.seed!(42)

N = 50_000
K1 = 1
K2 = 5
L = 4
J = 100

Y = rand(1:J,N)
X1 = randn(N,J,K1)
X2 = randn(J,K2)
D = randn(N,L)

α = randn(K1)
Π = randn(K2, L)
ξ_all = vcat(0, rand(J-1))
ξ = ξ_all[2:end]

Py = DiscreteChoice.y_prob(X1,X2,D, α, Π, ξ_all)
Y = zeros(Int64, N)
for i in eachindex(Y)
    Y[i] = sample(1:J, Weights(Py[i,:]))
end
[mean(Y.==i) for i in 1:J] |> minimum

DiscreteChoice.llk_sample(Y, X1, X2, D, α, Π, ξ)

#######################
# Estimate the model
#######################

model_homo = DCModel(Y,X1,X2,D; hetero_preference=false)
estDCModel!(model_homo; init_guess = vcat(α, vec(Π), ξ))
# estDCModel!(model_homo; init_guess = randn(K1+K2*L+J-1))
updateCoef!(model_homo)

@show model_homo.α .- α
@show model_homo.Π .- Π
@show model_homo.ξ .- ξ


text:
using Flux
using CuArrays

bn=BatchNorm(3,ϵ=Float32(1e-7),momentum=Float32(0.1)) |> gpu

@show typeof(bn.β)
@show typeof(bn.γ)
@show typeof(bn.λ)
@show typeof(bn.μ)
@show typeof(bn.σ)
@show typeof(bn.ϵ)
@show typeof(bn.momentum)

data = rand(28,28,3,2) |> gpu
bn(data)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
let
    # constructor
    term = Term("GO:0000001", "mitochondrion inheritance")
    @test isa(term, Term)
    @test term.id == 1
    @test term.name == "mitochondrion inheritance"
    @test !isobsolete(term)
    @test go"GO:0000001" == term
    @test go"0000001" == term
    @test string(term) == "Term(\"GO:0000001\", \"mitochondrion inheritance\")"

    obsolete_term = Term("GO:0000005", "ribosomal chaperone activity", obsolete=true)
    @test isobsolete(obsolete_term)
end


text:
struct ProblemParams{P,BT1,BT2}
    params :: P
    B₁ᵀ :: BT1
    B₂ :: BT2
end

function basic_unconstrained_problem(;tmax=1.0,iip=true)

  ns = 1
  y0 = 1.0
  α = 1.02
  p = α

  ode_rhs!(dy,y,p,t) = dy .= p*y
  ode_rhs(y,p,t) = p*y

  y₀ = ones(Float64,ns)
  u₀ = solvector(state=y₀)

  if iip
    f = ConstrainedODEFunction(ode_rhs!,_func_cache=u₀)
  else
    f = ConstrainedODEFunction(ode_rhs)
  end
  tspan = (0.0,1.0)
  prob = ODEProblem(f,u₀,tspan,p)

  xexact(t) = exp(α*t)

  return prob, xexact

end

function basic_unconstrained_if_problem(;tmax=1.0,iip=true)

  ns = 1
  y0 = 1.0
  α = 1.02
  p = α

  ode_rhs!(dy,y,p,t) = fill!(dy,0.0)
  ode_rhs(y,p,t) = zero(y)

  L = α*I

  y₀ = ones(Float64,ns)
  u₀ = solvector(state=y₀)

  if iip
    f = ConstrainedODEFunction(ode_rhs!,L,_func_cache=u₀)
  else
    f = ConstrainedODEFunction(ode_rhs,L)
  end
  tspan = (0.0,1.0)
  prob = ODEProblem(f,u₀,tspan,p)

  xexact(t) = exp(α*t)

  return prob, xexact

end

function basic_constrained_problem(;tmax=1.0,iip=true)

  U0 = 1.0
  g = 1.0
  α = 0.5
  ω = 5
  y₀ = Float64[0,0,U0,0]
  z₀ = Float64[0]

  params = [U0,g,α,ω];

  p₀ = ProblemParams(params,Array{Float64}(undef,4,1),Array{Float64}(undef,1,4));

  u₀ = solvector(state=y₀,constraint=z₀)
  du = deepcopy(u₀)

  function ode_rhs!(dy::Vector{Float64},y::Vector{Float64},p,t)
    dy .= 0.0
    dy[1] = y[3]
    dy[2] = y[4]
    dy[4] = -(y[1]-p.params[1]*t)*p.params[2]
    return dy
  end
  ode_rhs(y::Vector{Float64},p,t) = ode_rhs!(deepcopy(y₀),y,p,t)

  constraint_rhs!(dz::Vector{Float64},p,t) = dz .= Float64[p.params[1]]
  constraint_rhs(p,t) = constraint_rhs!(deepcopy(z₀),p,t)

  function op_constraint_force!(dy::Vector{Float64},z::Vector{Float64},p)
    @unpack B₁ᵀ = p
    dy .= B₁ᵀ*z
  end
  op_constraint_force(z::Vector{Float64},p) = op_constraint_force!(deepcopy(y₀),z,p)

  function constraint_op!(dz::Vector{Float64},y::Vector{Float64},p)
    @unpack B₂ = p
    dz .= B₂*y
  end
  constraint_op(y::Vector{Float64},p) = constraint_op!(deepcopy(z₀),y,p)

  function update_p!(q,u,p,t)
    y, z = state(u), constraint(u)
    @unpack B₁ᵀ, B₂ = q
    B₁ᵀ .= 0
    B₂ .= 0
    B₁ᵀ[3,1] = 1/(1+q.params[3]*sin(q.params[4]*t))
    B₂[1,3] = 1/(1+q.params[3]*sin(q.params[4]*t))
    return q
  end
  update_p(u,p,t) = update_p!(deepcopy(p),u,p,t)

  if iip
    f = ConstrainedODEFunction(ode_rhs!,constraint_rhs!,op_constraint_force!,
                              constraint_op!,_func_cache=deepcopy(du),
                                            param_update_func=update_p!)
  else
    f = ConstrainedODEFunction(ode_rhs,constraint_rhs,op_constraint_force,
                              constraint_op,param_update_func=update_p)
   end

  tspan = (0.0,tmax)
  p = deepcopy(p₀)
  prob = ODEProblem(f,u₀,tspan,p)

  yexact(t) = g*α*U0/ω*(-0.5*t^2 - cos(ω*t)/ω^2 + 1/ω^2)
  xexact(t) = U0*(t - α*cos(ω*t)/ω+α/ω)

  return prob, xexact, yexact
end

function cartesian_pendulum_problem(;tmax=1.0,iip=true)

  θ₀ = π/2
  l = 1.0
  g = 1.0
  y₀ = Float64[l*sin(θ₀),-l*cos(θ₀),0,0]
  z₀ = Float64[0.0, 0.0]

  u₀ = solvector(state=y₀,constraint=z₀)
  du = deepcopy(u₀)

  params = [l,g]
  p₀ = ProblemParams(params,Array{Float64}(undef,4,2),Array{Float64}(undef,2,4))

  function pendulum_rhs!(dy::Vector{Float64},y::Vector{Float64},p,t)
    dy .= 0.0
    dy[1] = y[3]
    dy[2] = y[4]
    dy[4] = -p.params[2]
    return dy
  end

  pendulum_rhs(y::Vector{Float64},p,t) = pendulum_rhs!(zero(y),y,p,t)

  length_constraint_rhs!(dz::Vector{Float64},p,t) = dz .= [0.0,p.params[1]^2]
  length_constraint_rhs(p,t) = length_constraint_rhs!(zero(z₀),p,t)

  function length_constraint_force!(dy::Vector{Float64},z::Vector{Float64},p)
    @unpack B₁ᵀ = p
    dy .= B₁ᵀ*z
  end
  length_constraint_force(z::Vector{Float64},p) = length_constraint_force!(zero(y₀),z,p)

  function length_constraint_op!(dz::Vector{Float64},y::Vector{Float64},p)
    @unpack B₂ = p
    dz .= B₂*y
  end
  length_constraint_op(y::Vector{Float64},p) = length_constraint_op!(zero(z₀),y,p)

  function update_p!(q,u,p,t)
    y, z = state(u), constraint(u)
    @unpack B₁ᵀ, B₂ = q
    fill!(B₁ᵀ,0.0)
    fill!(B₂,0.0)
    B₁ᵀ[3,1] = y[1]; B₁ᵀ[4,1] = y[2]; B₁ᵀ[1,2] = y[1]; B₁ᵀ[2,2] = y[2]
    B₂[1,3] = y[1]; B₂[1,4] = y[2]; B₂[2,1] = y[1]; B₂[2,2] = y[2]
    return q
  end
  update_p(u,p,t) = update_p!(deepcopy(p),u,p,t)

  if iip
    f = ConstrainedODEFunction(pendulum_rhs!,length_constraint_rhs!,length_constraint_force!,
                                length_constraint_op!,
                               _func_cache=deepcopy(du),param_update_func=update_p!)
  else
    f = ConstrainedODEFunction(pendulum_rhs,length_constraint_rhs,length_constraint_force,
                                length_constraint_op,param_update_func=update_p)
  end

  tspan = (0.0,tmax)
  p = deepcopy(p₀)
  prob = ODEProblem(f,u₀,tspan,p)

  # Get superconverged solution from the basic
  # problem expressed in theta
  function pendulum_theta(u,p,t)
      du = similar(u)
      du[1] = u[2]
      du[2] = -p^2*sin(u[1])
      du
  end

  u0 = [π/2,0.0]
  pex = g/l  # squared frequency
  tspan = (0.0,10.0)
  probex = ODEProblem(pendulum_theta,u0,tspan,pex)
  solex = solve(probex, Tsit5(), reltol=1e-16, abstol=1e-16);
  xexact(t) = sin(solex(t,idxs=1))
  yexact(t) = -cos(solex(t,idxs=1))

  return prob, xexact, yexact

end

function partitioned_problem(;tmax=1.0,iip=true)

  ω = 1.0
  βu = -0.2
  βv = -0.5

  par = [ω,βu,βv]

  U₀ = Float64[0,1]
  X₀ = Float64[1,0]
  Z₀ = Float64[0]
  u₀ = solvector(state=X₀,constraint=Z₀,aux_state=U₀)
  du = deepcopy(u₀)

  B₂ = Array{Float64}(undef,1,2)
  B₁ᵀ = Array{Float64}(undef,2,1)

  p₀ = ProblemParams(par,B₁ᵀ,B₂)

  L = Diagonal([βu,βv])

  function X_rhs!(dy,y,p,t)
    fill!(dy,0.0)
    return dy
  end
  X_rhs(y,p,t) = X_rhs!(deepcopy(X₀),y,p,t)

  function U_rhs!(dy,y,p,t)
    fill!(dy,0.0)
    ω = p.params[1]
    dy[1] = -ω*cos(ω*t)
    dy[2] = -ω*sin(ω*t)
    return dy
  end
  U_rhs(y,p,t) = U_rhs!(deepcopy(U₀),y,p,t)

  ode_rhs! = ArrayPartition((X_rhs!,U_rhs!))
  ode_rhs = ArrayPartition((X_rhs,U_rhs))

  constraint_rhs!(dz,p,t) = dz .= Float64[0]
  constraint_rhs(p,t) = constraint_rhs!(deepcopy(Z₀),p,t)

  function op_constraint_force!(dy,z,p)
    @unpack B₁ᵀ = p
    dy .= B₁ᵀ*z
    return dy
  end
  op_constraint_force(z,p) = op_constraint_force!(deepcopy(X₀),z,p)

  function constraint_op!(dz,y,p)
    @unpack B₂ = p
    dz .= B₂*y
  end
  constraint_op(y,p) = constraint_op!(deepcopy(Z₀),y,p)

  function update_p!(q,u,p,t)
    x = aux_state(u)
    @unpack B₁ᵀ, B₂ = q

    B₁ᵀ[1,1] = x[1]
    B₁ᵀ[2,1] = x[2]
    B₂[1,1] = x[1]
    B₂[1,2] = x[2]
    return q
  end
  update_p(u,p,t) = update_p!(deepcopy(p),u,p,t)

  if iip
    f = ConstrainedODEFunction(ode_rhs!,constraint_rhs!,op_constraint_force!,
                              constraint_op!,L,_func_cache=deepcopy(du),
                              param_update_func=update_p!)
  else
    f = ConstrainedODEFunction(ode_rhs,constraint_rhs,op_constraint_force,
                              constraint_op,L,param_update_func=update_p)
  end

  tspan = (0.0,tmax)
  p = deepcopy(p₀)
  update_p!(p,u₀,p,0.0)
  prob = ODEProblem(f,u₀,tspan,p)

  # function fex(du,u,p,t)
  #   UV = u.x[1]
  #   xy = u.x[2]
  #   ω = p[1]
  #   βu = p[2]
  #   βv = p[3]
  #   du[1] = -ω*cos(ω*t)
  #   du[2] = -ω*sin(ω*t)
  #   Usq = u[1]^2+u[2]^2
  #   du[3] = (βu-u[1]/Usq*(du[1]+βu*u[1]))*u[3] - u[1]/Usq*(du[2]+βv*u[2])*u[4]
  #   du[4] = -u[2]/Usq*(du[1]+βu*u[1])*u[3] + (βv-u[2]/Usq*(du[2]+βv*u[2]))*u[4]
  #   return nothing
  # end
  #
  # tspan = (0.0,tmax)
  # u₀ex = SaddleVector(U₀,X₀)
  #
  # probex = ODEProblem(fex,u₀ex,tspan,par)
  # solex = solve(probex, Tsit5(), reltol=1e-16, abstol=1e-16)
  # xexact(t) = solex(t,idxs=3)
  # yexact(t) = solex(t,idxs=4)

  fex(t) = exp(0.5*(βu+βv)*t)*exp(0.25/ω*sin(2ω*t)*(βu-βv))
  xexact(t) = cos(ω*t)*fex(t)
  yexact(t) = sin(ω*t)*fex(t)

  return prob, xexact, yexact

end


text:
# compilation backend that uses CUDD

using CUDD

mutable struct CuddMgr <: DiceManager
    cuddmgr::Ptr{Nothing}
    probs::Dict{Int,Float64}
end

function CuddMgr()
    cudd_mgr = initialize_cudd()
    Cudd_DisableGarbageCollection(cudd_mgr) # note: still need to ref because CUDD can delete nodes without doing a GC pass
    mgr = CuddMgr(cudd_mgr, Dict{Int,Float64}())
    finalizer(mgr) do x
        Cudd_Quit(x.cuddmgr)
    end
end

##################################
# core functionality
##################################

constant(mgr::CuddMgr, c:: Bool) =
    c ? Cudd_ReadOne(mgr.cuddmgr) : Cudd_ReadLogicZero(mgr.cuddmgr)

biconditional(mgr::CuddMgr, x, y) =
    rref(Cudd_bddXnor(mgr.cuddmgr, x, y))

conjoin(mgr::CuddMgr, x, y) =
    rref(Cudd_bddAnd(mgr.cuddmgr, x, y))

disjoin(mgr::CuddMgr, x, y) =
    rref(Cudd_bddOr(mgr.cuddmgr, x, y))

negate(::CuddMgr, x) =
    Cudd_Not(x)

ite(mgr::CuddMgr, cond, then, elze) =
    rref(Cudd_bddIte(mgr.cuddmgr, cond, then, elze))

new_var(mgr::CuddMgr, prob) = begin
    x = rref(Cudd_bddNewVar(mgr.cuddmgr))
    mgr.probs[decisionvar(mgr, x)] = prob
    x
end

function infer(mgr::CuddMgr, x)

    cache = Dict{Tuple{Ptr{Nothing},Bool},Float64}()
    t = constant(mgr, true)
    cache[(t,false)] = log(one(Float64))
    cache[(t,true)] = log(zero(Float64))

    rec(y, c) =
        if Cudd_IsComplement(y)
            rec(Cudd_Regular(y), !c)
        else get!(cache, (y,c)) do
                v = decisionvar(mgr,y)
                prob = mgr.probs[v]
                a = log(prob) + rec(Cudd_T(y), c)
                b = log(1.0-prob) + rec(Cudd_E(y), c)
                if (!isfinite(a))
                    b
                elseif (!isfinite(b))
                    a
                else
                    max(a,b) + log1p(exp(-abs(a-b)))
                end
            end
        end

    logprob = rec(x, false)
    exp(logprob)
end

##################################
# additional CUDD-based functionality
##################################

function Base.show(io::IO, mgr::CuddMgr, x)
    if !issat(mgr, x)
        print(io, "(false)")
    elseif isvalid(mgr, x)
        print(io, "(true)")
    elseif isposliteral(mgr, x)
        print(io, "(f$(decisionvar(mgr, x)))")
    elseif isnegliteral(mgr, x)
        print(io, "(-f$(decisionvar(mgr, x)))")
    else
        print(io, "@$(hash(x)÷ 10000000000000)")
    end
end

function Base.show(io::IO, x::CuddMgr)
    print(io, "$(typeof(x))@$(hash(x)÷ 10000000000000)")
end

isconstant(x) =
    isone(Cudd_IsConstant(x))

isliteral(x::DistBool) =
    isliteral(x.mgr, x.bit)

isliteral(::CuddMgr, x) =
    (!isconstant(x) &&
     isconstant(Cudd_T(x)) &&
     isconstant(Cudd_E(x)))

isposliteral(x::DistBool) =
    isposliteral(x.mgr, x.bit)

isposliteral(mgr::CuddMgr, x) =
    isliteral(mgr,x) &&
    (x === Cudd_bddIthVar(mgr.cuddmgr, decisionvar(mgr,x)))

isnegliteral(x::DistBool) =
    isnegliteral(x.mgr, x.bit)

isnegliteral(mgr::CuddMgr, x) =
    isliteral(mgr,x) &&
    (x !== Cudd_bddIthVar(mgr.cuddmgr, decisionvar(mgr,x)))

issat(x::DistBool) =
    issat(x.mgr, x.bit)

issat(mgr::CuddMgr, x) =
    x !== constant(mgr, false)

isvalid(x::DistBool) =
    isvalid(x.mgr, x.bit)

isvalid(mgr::CuddMgr, x) =
    x === constant(mgr, true)

num_nodes(bits::Vector{DistBool}; as_add=true) =
    num_nodes(bits[1].mgr, map(b -> b.bit, bits); as_add)

num_nodes(x; as_add=true) =
    num_nodes(bools(x); as_add)

num_nodes(mgr::CuddMgr, xs::Vector{<:Ptr}; as_add=true) = begin
    as_add && (xs = map(x -> rref(Cudd_BddToAdd(mgr.cuddmgr, x)), xs))
    Cudd_SharingSize(xs, length(xs))
end

num_flips(bits::Vector{DistBool}) =
    num_vars(bits[1].mgr, map(b -> b.bit, bits))

num_flips(x) =
    num_flips(bools(x))

num_vars(mgr::CuddMgr, xs::Vector{<:Ptr}) = begin
    Cudd_VectorSupportSize(mgr.cuddmgr, xs, length(xs))
end

num_vars(mgr::CuddMgr) =
    Cudd_ReadSize(mgr.cuddmgr)

decisionvar(_::CuddMgr, x) =
    Cudd_NodeReadIndex(x)

dump_dot(bits::Vector{DistBool}, filename; as_add=true) =
    dump_dot(bits[1].mgr, map(b -> b.bit, bits), filename; as_add)

dump_dot(x, filename; as_add=true) =
    dump_dot(bools(x), filename; as_add)

mutable struct FILE end

dump_dot(mgr::CuddMgr, xs::Vector{<:Ptr}, filename; as_add=true) = begin
    # convert to ADDs in order to properly print terminals
    if as_add
        xs = map(x -> rref(Cudd_BddToAdd(mgr.cuddmgr, x)), xs)
    end
    outfile = ccall(:fopen, Ptr{FILE}, (Cstring, Cstring), filename, "w")
    Cudd_DumpDot(mgr.cuddmgr, length(xs), xs, C_NULL, C_NULL, outfile)
    @assert ccall(:fclose, Cint, (Ptr{FILE},), outfile) == 0
    nothing
end

##################################
# CUDD Utilities
##################################

rref(x) = begin
    ref(x)
    x
end


text:
module FEMMCorotBeamModule

using LinearAlgebra: norm, Transpose, mul!
using FinEtools
using FinEtools.IntegDomainModule: IntegDomain
import FinEtoolsDeforLinear.MatDeforElastIsoModule: MatDeforElastIso
using ..FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!, MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA

"""
    FEMMDeforLinear{MR<:AbstractDeforModelRed,  S<:AbstractFESet, F<:Function, M<:AbstractMatDeforLinearElastic} <: FEMMCorotBeam

Class for linear deformation finite element modeling machine.
"""
mutable struct FEMMCorotBeam{S<:AbstractFESet, F<:Function} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # integration domain data
    material::MatDeforElastIso # material object
    # The attributes below are buffers used in various operations.
    _ecoords0::FFltMat
    _ecoords1::FFltMat
    _edisp1::FFltMat
    _evel1::FFltMat
    _evel1f::FFltMat
    _dofnums::FIntMat
    _F0::FFltMat
    _Ft::FFltMat
    _FtI::FFltMat
    _FtJ::FFltMat
    _Te::FFltMat
    _tempelmat1::FFltMat
    _tempelmat2::FFltMat
    _tempelmat3::FFltMat
    _elmat::FFltMat
    _elmatTe::FFltMat
    _elmato::FFltMat
    _elvec::FFltVec
    _elvecf::FFltVec
    _aN::FFltMat
    _dN::FFltVec
    _DN::FFltMat
    _PN::FFltVec
    _LF::FFltVec
    _RI::FFltMat
    _RJ::FFltMat
    _OS::FFltMat
end

function FEMMCorotBeam(integdomain::IntegDomain{S, F}, material::MatDeforElastIso) where {S<:FESetL2CorotBeam, F<:Function}
    _ecoords0 = fill(0.0, 2, 3);
    _ecoords1 = fill(0.0, 2, 3)
    _edisp1 = fill(0.0, 2, 3);
    _evel1 = fill(0.0, 2, 6);
    _evel1f = fill(0.0, 2, 6)
    _dofnums = zeros(FInt, 1, 12);
    _F0 = fill(0.0, 3, 3);
    _Ft = fill(0.0, 3, 3);
    _FtI = fill(0.0, 3, 3);
    _FtJ = fill(0.0, 3, 3)
    _Te = fill(0.0, 12, 12)
    _tempelmat1 = fill(0.0, 12, 12);
    _tempelmat2 = fill(0.0, 12, 12);
    _tempelmat3 = fill(0.0, 12, 12)
    _elmat = fill(0.0, 12, 12);
    _elmatTe = fill(0.0, 12, 12);
    _elmato = fill(0.0, 12, 12)
    _elvec = fill(0.0, 12);
    _elvecf = fill(0.0, 12)
    _aN = fill(0.0, 6, 12)
    _dN = fill(0.0, 6)
    _DN = fill(0.0, 6, 6)
    _PN = fill(0.0, 6)
    _LF = fill(0.0, 12)
    _RI = fill(0.0, 3, 3);
    _RJ = fill(0.0, 3, 3);
    _OS = fill(0.0, 3, 3)
    return FEMMCorotBeam(integdomain, material,
     _ecoords0, _ecoords1, _edisp1, _evel1, _evel1f,
     _dofnums,
     _F0, _Ft, _FtI, _FtJ, _Te,
     _tempelmat1, _tempelmat2, _tempelmat3, _elmat, _elmatTe, _elmato,
     _elvec, _elvecf,
     _aN, _dN, _DN, _PN, _LF,
     _RI, _RJ, _OS)
end

function _transfmat!(Te, Ft)
    Te[1:3, 1:3] = Te[4:6, 4:6] = Te[7:9, 7:9] = Te[10:12, 10:12] = Ft
    return Te
end

"""
    mass(self::FEMMCorotBeam,  assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the consistent mass matrix

This is a general routine for the abstract linear-deformation  FEMM.
"""
function mass(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}; mass_type=MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    ecoords0, ecoords1, edisp1, dofnums = self._ecoords0, self._ecoords1, self._edisp1, self._dofnums
    F0, Ft, FtI, FtJ, Te = self._F0, self._Ft, self._FtI, self._FtJ, self._Te
    R1I, R1J = self._RI, self._RJ
    elmat, elmatTe = self._elmat, self._elmatTe
    dN = self._dN
    rho = massdensity(self.material)
    A, I1, I2, I3, x1x2_vector = fes.A, fes.I1, fes.I2, fes.I3, fes.x1x2_vector
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords0, fes.conn[i]);
        gathervalues_asmat!(u1, edisp1, fes.conn[i]);
        ecoords1 .= ecoords0 .+ edisp1
        R1I[:] .= Rfield1.values[fes.conn[i][1], :];
        R1J[:] .= Rfield1.values[fes.conn[i][2], :];
        fill!(elmat,  0.0); # Initialize element matrix
        L1, Ft, dN = local_frame_and_def!(Ft, dN, F0, FtI, FtJ, ecoords0, x1x2_vector[i], ecoords1, R1I, R1J);
        _transfmat!(Te, Ft)
        L0 = norm(ecoords0[2,:]-ecoords0[1,:]);
        local_mass!(elmat, A[i], I1[i], I2[i], I3[i], rho, L0, mass_type);
        mul!(elmatTe, elmat, Transpose(Te))
        mul!(elmat, Te, elmatTe)
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function mass(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}; mass_type=MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return mass(self, assembler, geom0, u1, Rfield1, dchi; mass_type = mass_type);
end

"""
    gyroscopic(self::FEMMCorotBeam,  assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the quadratic-inertial-term (gyroscopic) mass matrix

This is a general routine for the abstract linear-deformation  FEMM.
"""
function gyroscopic(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, v1::NodalField{T}, dchi::NodalField{TI}; mass_type=MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    ecoords0, ecoords1, edisp1, dofnums = self._ecoords0, self._ecoords1, self._edisp1, self._dofnums
    F0, Ft, FtI, FtJ, Te = self._F0, self._Ft, self._FtI, self._FtJ, self._Te
    R1I, R1J = self._RI, self._RJ
    elmat, elmatTe = self._elmat, self._elmatTe
    dN = self._dN
    evel1, evel1f = self._evel1, self._evel1f
    Ge1, Ge2, Ge = self._tempelmat1, self._tempelmat2, self._tempelmat3
    OmegaTilde = self._elmato
    OS = self._OS
    rho = massdensity(self.material)
    A, I1, I2, I3, x1x2_vector = fes.A, fes.I1, fes.I2, fes.I3, fes.x1x2_vector
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords0, fes.conn[i]);
        gathervalues_asmat!(u1, edisp1, fes.conn[i]);
        ecoords1 .= ecoords0 .+ edisp1
        R1I[:] .= Rfield1.values[fes.conn[i][1], :];
        R1J[:] .= Rfield1.values[fes.conn[i][2], :];
        fill!(elmat,  0.0); # Initialize element matrix
        L1, Ft, dN = local_frame_and_def!(Ft, dN, F0, FtI, FtJ, ecoords0, x1x2_vector[i], ecoords1, R1I, R1J);
        _transfmat!(Te, Ft)
        L0 = norm(ecoords0[2,:]-ecoords0[1,:]);
        local_mass!(elmat, A[i], I1[i], I2[i], I3[i], rho, L0, mass_type);
        mul!(elmatTe, elmat, Transpose(Te))
        mul!(elmat, Te, elmatTe)
        gathervalues_asmat!(v1, evel1, fes.conn[i]);
        evel1f[:, 1:3] = evel1[:,1:3]*Ft
        evel1f[:, 4:6] = evel1[:,4:6]*Ft
        Omega = (evel1f[1,4]+evel1f[2,4])/2*Ft[:,1] + (evel1f[1,3]-evel1f[2,3])/L1*Ft[:,2] + (evel1f[2,2]-evel1f[1,2])/L1*Ft[:,3];
        skewmat!(OS, Omega);
        _transfmat!(OmegaTilde, OS)
        mul!(Ge1, OmegaTilde, elmat)
        mul!(Ge2, elmat, OmegaTilde)
        @. Ge = Ge1 - Ge2;
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, Ge, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function gyroscopic(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, v1::NodalField{T}, dchi::NodalField{TI}; mass_type=MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return gyroscopic(self, assembler, geom0, u1, Rfield1, v1, dchi; mass_type = mass_type);
end

"""
    stiffness(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{T}) where {ASS<:AbstractSysmatAssembler, T<:Number}

Compute the material stiffness matrix.
"""
function stiffness(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    ecoords0, ecoords1, edisp1, dofnums = self._ecoords0, self._ecoords1, self._edisp1, self._dofnums
    F0, Ft, FtI, FtJ, Te = self._F0, self._Ft, self._FtI, self._FtJ, self._Te
    R1I, R1J = self._RI, self._RJ
    elmat, elmatTe = self._elmat, self._elmatTe
    aN, dN, DN = self._aN, self._dN, self._DN
    E = self.material.E
    G = E / 2 / (1 + self.material.nu)::Float64
    A, I2, I3, J, A2s, A3s, x1x2_vector = fes.A, fes.I2, fes.I3, fes.J, fes.A2s, fes.A3s, fes.x1x2_vector
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords0, fes.conn[i]);
        gathervalues_asmat!(u1, edisp1, fes.conn[i]);
        ecoords1 .= ecoords0 .+ edisp1
        R1I[:] .= Rfield1.values[fes.conn[i][1], :];
        R1J[:] .= Rfield1.values[fes.conn[i][2], :];
        fill!(elmat,  0.0); # Initialize element matrix
        L1, Ft, dN = local_frame_and_def!(Ft, dN, F0, FtI, FtJ, ecoords0, x1x2_vector[i], ecoords1, R1I, R1J);
        _transfmat!(Te, Ft)
        local_stiffness!(elmat, E, G, A[i], I2[i], I3[i], J[i], A2s[i], A3s[i], L1, aN, DN);
        mul!(elmatTe, elmat, Transpose(Te))
        mul!(elmat, Te, elmatTe)
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function stiffness(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom0, u1, Rfield1, dchi);
end

"""
    geostiffness(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{T}) where {ASS<:AbstractSysmatAssembler, T<:Number}

Compute the geometric stiffness matrix.
"""
function geostiffness(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    ecoords0, ecoords1, edisp1, dofnums = self._ecoords0, self._ecoords1, self._edisp1, self._dofnums
    F0, Ft, FtI, FtJ, Te = self._F0, self._Ft, self._FtI, self._FtJ, self._Te
    R1I, R1J = self._RI, self._RJ
    elmat, elmatTe = self._elmat, self._elmatTe
    aN, dN, DN, PN = self._aN, self._dN, self._DN, self._PN
    E = self.material.E
    G = E / 2 / (1 + self.material.nu)
    A, I2, I3, J, A2s, A3s, x1x2_vector = fes.A, fes.I2, fes.I3, fes.J, fes.A2s, fes.A3s, fes.x1x2_vector
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords0, fes.conn[i]);
        gathervalues_asmat!(u1, edisp1, fes.conn[i]);
        ecoords1 .= ecoords0 .+ edisp1
        R1I[:] .= Rfield1.values[fes.conn[i][1], :];
        R1J[:] .= Rfield1.values[fes.conn[i][2], :];
        fill!(elmat,  0.0); # Initialize element matrix
        L1, Ft, dN = local_frame_and_def!(Ft, dN, F0, FtI, FtJ, ecoords0, x1x2_vector[i], ecoords1, R1I, R1J);
        _transfmat!(Te, Ft)
        natural_forces!(PN, E, G, A[i], I2[i], I3[i], J[i], A2s[i], A3s[i], L1, dN, DN)
        local_geometric_stiffness!(elmat, A[i], I2[i], I3[i], PN, L1)
        mul!(elmatTe, elmat, Transpose(Te))
        mul!(elmat, Te, elmatTe)
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function geostiffness(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return geostiffness(self, assembler, geom0, u1, Rfield1, dchi);
end

"""
    restoringforce(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{T}) where {ASS<:AbstractSysmatAssembler, T<:Number}

Compute the vector of the restoring elastic forces
"""
function restoringforce(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {ASS<:AbstractSysvecAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    ecoords0, ecoords1, edisp1, dofnums = self._ecoords0, self._ecoords1, self._edisp1, self._dofnums
    F0, Ft, FtI, FtJ, Te = self._F0, self._Ft, self._FtI, self._FtJ, self._Te
    R1I, R1J = self._RI, self._RJ
    elmat, elmatTe = self._elmat, self._elmatTe
    aN, dN, DN, PN, LF = self._aN, self._dN, self._DN, self._PN, self._LF
    elvec = self._elvec
    E = self.material.E
    G = E / 2 / (1 + self.material.nu)
    A, I2, I3, J, A2s, A3s, x1x2_vector = fes.A, fes.I2, fes.I3, fes.J, fes.A2s, fes.A3s, fes.x1x2_vector
    startassembly!(assembler, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords0, fes.conn[i]);
        gathervalues_asmat!(u1, edisp1, fes.conn[i]);
        ecoords1 .= ecoords0 .+ edisp1
        R1I[:] .= Rfield1.values[fes.conn[i][1], :];
        R1J[:] .= Rfield1.values[fes.conn[i][2], :];
        fill!(elmat,  0.0); # Initialize element matrix
        L1, Ft, dN = local_frame_and_def!(Ft, dN, F0, FtI, FtJ, ecoords0, x1x2_vector[i], ecoords1, R1I, R1J);
        _transfmat!(Te, Ft)
        natural_forces!(PN, E, G, A[i], I2[i], I3[i], J[i], A2s[i], A3s[i], L1, dN, DN)
        local_forces!(LF, PN, L1, aN)
        mul!(elvec, Te, -LF)
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, elvec, dofnums);
    end # Loop over elements
    return makevector!(assembler);
end

function restoringforce(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {T<:Number, TI<:Number}
    assembler = SysvecAssembler();
    return restoringforce(self, assembler, geom0, u1, Rfield1, dchi);
end

"""
    distribloads_global(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{T}, fi) where {T<:Number}

Compute the load vector due to distributed loads.

Compute the global load vector corresponding to applied distributed
load. Here it means force per unit length of the beam,
in the configuration u1,Rfield1.

Note: the force intensity must be uniform across the entire element.
Note: the force intensity is given in the global coordinates.
"""
function distribloads_global(self::FEMMCorotBeam, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}, fi) where {ASS<:AbstractSysvecAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    ecoords0, ecoords1, edisp1, dofnums = self._ecoords0, self._ecoords1, self._edisp1, self._dofnums
    F0, Ft, FtI, FtJ, Te = self._F0, self._Ft, self._FtI, self._FtJ, self._Te
    R1I, R1J = self._RI, self._RJ
    elmat, elmatTe = self._elmat, self._elmatTe
    aN, dN, DN, PN = self._aN, self._dN, self._DN, self._PN
    elvec, elvecf = self._elvec, self._elvecf
    Lforce = fill(0.0, 3)
    ignore = fill(0.0 , 0, 0)
    E = self.material.E
    G = E / 2 / (1 + self.material.nu)
    A, I2, I3, J, x1x2_vector = fes.A, fes.I2, fes.I3, fes.J, fes.x1x2_vector
    startassembly!(assembler, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords0, fes.conn[i]);
        gathervalues_asmat!(u1, edisp1, fes.conn[i]);
        ecoords1 .= ecoords0 .+ edisp1
        R1I[:] .= Rfield1.values[fes.conn[i][1], :];
        R1J[:] .= Rfield1.values[fes.conn[i][2], :];
        fill!(elmat,  0.0); # Initialize element matrix
        L1, Ft, dN = local_frame_and_def!(Ft, dN, F0, FtI, FtJ, ecoords0, x1x2_vector[i], ecoords1, R1I, R1J);
        _transfmat!(Te, Ft)
        L0 = norm(ecoords0[2,:]-ecoords0[1,:]);
        force = updateforce!(fi, ignore, ignore, fes.label[i]); # retrieve the applied load
        Lforce = Ft' * force # local distributed load components
        elvecf[1] = Lforce[1]*L0/2;
        elvecf[2] = Lforce[2]*L0/2;
        elvecf[3] = Lforce[3]*L0/2;
        elvecf[4] = 0;
        elvecf[5] = -Lforce[3]*L0^2/12;
        elvecf[6] = +Lforce[2]*L0^2/12;
        elvecf[7] = Lforce[1]*L0/2;
        elvecf[8] = Lforce[2]*L0/2;
        elvecf[9] = Lforce[3]*L0/2;
        elvecf[10] = 0;
        elvecf[11] = +Lforce[3]*L0^2/12;
        elvecf[12] = -Lforce[2]*L0^2/12;
        mul!(elvec, Te, elvecf)
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, elvec, dofnums);
    end # Loop over elements
    return makevector!(assembler);
end

function distribloads_global(self::FEMMCorotBeam, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}, fi) where {T<:Number, TI<:Number}
    assembler = SysvecAssembler();
    return distribloads_global(self, assembler, geom0, u1, Rfield1, dchi, fi);
end

end # module


text:
#=
Each character on a computer is assigned a unique code and the preferred standard is ASCII (American Standard Code for Information Interchange). For example, uppercase A = 65, asterisk (*) = 42, and lowercase k = 107.

A modern encryption method is to take a text file, convert the bytes to ASCII, then XOR each byte with a given value, taken from a secret key. The advantage with the XOR function is that using the same encryption key on the cipher text, restores the plain text; for example, 65 XOR 42 = 107, then 107 XOR 42 = 65.

For unbreakable encryption, the key is the same length as the plain text message, and the key is made up of random bytes. The user would keep the encrypted message and the encryption key in different locations, and without both "halves", it is impossible to decrypt the message.

Unfortunately, this method is impractical for most users, so the modified method is to use a password as a key. If the password is shorter than the message, which is likely, the key is repeated cyclically throughout the message. The balance for this method is using a sufficiently long password key for security, but short enough to be memorable.

Your task has been made easy, as the encryption key consists of three lower case characters. Using cipher.txt (right click and 'Save Link/Target As...'), a file containing the encrypted ASCII codes, and the knowledge that the plain text must contain common English words, decrypt the message and find the sum of the ASCII values in the original text.
=#
using DelimitedFiles

function check(cypher, key)
  y = UInt8[]
  for i in 1:length(cypher)
    push!(y, xor(cypher[i], key[mod(i,length(key))+1]))
  end
  result = String(y)

  # https://en.wikipedia.org/wiki/Most_common_words_in_English
  commonwords = split("the|be|to|of|and|a|in|that|have|i|it|for|not|on|with|he|as|you|do|at|this|but|his|by|from|they|we|say|her|she|or|an|will|my|one", '|')
  s = filter((x)->occursin(x, result), commonwords)
  if length(s) > 20
    return sum(codeunits(result))
  end
  return 0
end

function calc()
  encryptedtext = vec(readdlm(dirname(@__FILE__()) * "/../Resources/p59.txt", ',', UInt8))
  for a in 'a':'z', b in 'a':'z', c in 'a':'z'
    key = codeunits("$a$b$c")
    r = check(encryptedtext, key)
    if r > 0
      return r
    end
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using StatsBase
using Plots.PlotMeasures

include("models/data.jl")
include("models/data_kw.jl")

d=mordat

function plot_sterbenszeiten(mordat; kw...)
    rs = unique(mordat.date)
    plot_unstacked_heatmap(mordat;
                           labelf = x -> let (g,a) = x; (g!="Männlich" ? (alterstring(a)*"    ") : "") * label_geschlecht[g]; end,
                           ylabel = "Alter   ♂♀",
                           xticks = (1:52:length(rs), string.(2000:1:2021)),
                           xrotation = -90,
                           xmirror=true,
                           c = tot_cgrad,
                           top_margin = 20px,
                           bottom_margin = 20px,
                           kw...
                               ) do p,d
                                   p
                               end
end

# @transform! groupby(mordat, [:alter, :geschlecht]) :PD_12months = gliding_mean(:PD),

@transform!(groupby(mordat, [:alter, :geschlecht]),
            :PD_rel2000 = :PD ./ first(skipmissing(:PD_12months)),
            :PD_relmean = :PD ./ mean(skipmissing(:PD)),
#            :PD_rel_12months = :PD_12months ./ first(skipmissing(:PD_12months)),
#            :PD_relmean_12months = :PD_12months ./ mean(skipmissing(:PD_12months)),
            )

@transform!(mordat,
            :row = collect(zip(:geschlecht, :alter)),
            :col = :date,
            :value = 100 .* :PD_relmean
            )

# filter(x -> x.date>Date(2016,1,1), mordat)
p = plot_sterbenszeiten(
    mordat,
    title = "Rythmen der Übersterblichkeit in %\nBundesrepublik Deutschland"
)
DeutschlandDaten.savefigs(p, "relative_Übersterblichkeit")

p = plot_sterbenszeiten(
    mordat,
    title = "Übersterblichkeit der Geburtsjahrgänge in %\nBundesrepublik Deutschland")
# Kohorten-Linien
kohdat = filter(x-> rem(x.geburt,10)==0 ,
                DataFrame([ (geburt = j - x.start, woche = 1+52*(j-2000), alter_von = i*2-1.5)
                            for (i,x) in enumerate(DeutschlandDaten.Altersgruppen_KW)
                                for j in 2000:2021 ]))
@df kohdat plot!(:woche, :alter_von, group=:geburt,
                 c=:white,
                 legend=nothing)#, bg_inside = nothing,
                 # xlim=nothing,
                 #ylim=nothing,
                 #framestyle=nothing)
annotate!([ ( x[end,:woche], x[end,:alter_von], ("* $(x[end,:geburt])", 7, :left, :bottom, :white)) for x in groupby(kohdat, :geburt)
               if size(x)[1]>1])
#plotquelle!(p; subplot = 4) # breaks plot!

DeutschlandDaten.savefigs(p,"relative_Übersterblichkeit_kohorten")

@transform!(mordat,
            :col = collect(zip(:jahr, :geschlecht)),
            :row = :kw
            )

for alter=DeutschlandDaten.Altersgruppen_KW
    title = "Sterberate Deutschland im Alter $alter von 10'000"
    d = filter(x -> x.alter==alter, mordat)
    us = unstack(d[:,[:row,:col,:PD]], :col, :PD)
    p = heatmap((Matrix(us[:,2:end]));
                title = title,
                yticks = 1:5:52,
                xticks = (vcat(1:5:22, 23:5:44), repeat(collect(string.(2000:5:2021)),2)),
                xrotation = -90,
                c = tot_cgrad )
    plot!([22.5]; seriestype=:vline, color=:black, legend=nothing, width=2)
    annotate!([(12,45,"Männlich"), (23+12,45,"Weiblich")])
    DeutschlandDaten.savefigs(p, "Sterberate_$(alter)")
end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
using CUDAdrv
using Base.Test

dev = CuDevice(0)
ctx = CuContext(dev)

md = CuModuleFile(joinpath(dirname(@__FILE__), "vadd.ptx"))
vadd = CuFunction(md, "kernel_vadd")

dims = (3,4)
a = round(rand(Float32, dims) * 100)
b = round(rand(Float32, dims) * 100)

d_a = CuArray(a)
d_b = CuArray(b)
d_c = CuArray(Float32, dims)

len = prod(dims)
cudacall(vadd, len, 1, (DevicePtr{Cfloat},DevicePtr{Cfloat},DevicePtr{Cfloat}), d_a, d_b, d_c)
c = Array(d_c)
@test a+b ≈ c

destroy(ctx)


text:
using Genie.Router
using Genie.Requests
using ItemsController: items

route("/") do
  serve_static_file("welcome.html")
end

route("/items", () -> items(Val(:view), Val(:GET)); method = GET, named = :get_items)
route("/items", () -> items(Val(:view), Val(:POST)); method = POST, named = :post_items)

route("/items/:id", () -> items(Val(:view), Val(:GET), payload(:id)); method = GET, named = :get_items_id)
route("/items/:id", () -> items(Val(:view), Val(:POST), payload(:id)); method = POST, named = :post_items_id)

route("/api/items", () -> items(Val(:api), Val(:GET)); method = GET, named = :get_api_items)
route("/api/items", () -> items(Val(:api), Val(:POST)); method = POST, named = :post_api_items)

route("/api/items/:id", () -> items(Val(:api), Val(:GET), payload(:id)); method = GET, named = :get_api_items_id)
route("/api/items/:id", () -> items(Val(:api), Val(:PUT), payload(:id)); method = PUT, named = :put_api_items_id)
route("/api/items/:id", () -> items(Val(:api), Val(:DELETE), payload(:id)); method = DELETE, named = :delete_api_items_id)

max_stars_repo_name: jaantollander/genie-webapp-db

max_issues_repo_name: jaantollander/genie-webapp-db

max_forks_repo_name: jaantollander/genie-webapp-db


text:
# Uncomment this when all CUDA modules support precompilation
# isdefined(Base, :__precompile__) && __precompile__()
#module Knet
using Compat
info("Loading Knet")

### GPU detection and initialization
include("util/gpu.jl");		export gpu, @gpu, @useifgpu, setseed
@useifgpu CUDArt
@useifgpu CUBLAS
@useifgpu CUSPARSE
@useifgpu CUDNN

### GPU dependent includes
@gpu include("util/cudart.jl");		# export copysync!, fillsync!
@gpu include("util/curand.jl");
@gpu include("util/cusparse.jl");
@gpu include("util/deepcopy.jl"); 	# export cpucopy, gpucopy
@gpu include("util/linalg.jl");

### GPU independent utilities
include("util/dbg.jl");		export @date, @dbg, gpusync
include("util/rgen.jl");	export Gaussian, Uniform, Constant, Identity, Xavier, Bernoulli, GlorotUniform, GlorotNormal
include("util/array.jl");	# export isapprox
include("util/colops.jl");	export csize, clength, ccount, csub, cget, size2, minibatch
include("util/conv_pool_cpu.jl");

### Main
include("op.jl")
include("op/actf.jl")
include("op/broadcast.jl")
include("op/conv.jl")
include("op/dot.jl")
include("op/lrn.jl")
include("op/nce.jl")
include("op/pool.jl")
include("op/genop.jl")
include("net.jl");              export setp, wnorm, gnorm, reset!, clean, stack_isempty #export Reg, Net, set!, inc!, registers, params, ninputs, out, dif, stack_empty!,
include("net/forw.jl");         export forw, sforw
include("net/back.jl");         export back, sback
include("net/initforw.jl")
include("net/initback.jl")
include("compiler.jl");		export @knet, compile, _comp_parse_def # @knet needs _comp_parse_def
include("gradcheck.jl");  	export gradcheck
include("update.jl");		export update!
include("loss.jl");		export quadloss, softloss, zeroone, xentloss # TODO-TEST: logploss, xentloss, percloss, scalloss,
#include("kfun.jl")

# To be deprecated:
include("data/ItemTensor.jl");		export ItemTensor
include("data/S2SData.jl");     	export S2SData, maxtoken
include("data/SequencePerLine.jl"); 	export SequencePerLine
include("data/SketchEngine.jl"); 	export SketchEngine
include("data/TagData.jl"); 		export TagData, sponehot

# Knet8 stuff
include("tmplike.jl");		export tmplike, tmpdict, tmpkeep

# Load kernels from CUDArt
function __init__()
# Let's just use one gpu for now.
#    @gpu CUDArt.init(CUDArt.devices(i->true))
    @gpu CUDArt.init(CUDArt.device())
end

__init__()
info("Loading Knet done.")
#end # module


text:
#####
##### Tracer advection operator
#####

"""
    div_uc(i, j, k, grid, advection, U, c)

Calculates the divergence of the flux of a tracer quantity c being advected by
a velocity field U = (u, v, w), ∇·(Uc),

    1/V * [δxᶜᵃᵃ(Ax * u * ℑxᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * v * ℑyᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * w * ℑzᵃᵃᶠ(c))]

which will end up at the location `ccc`.
"""
@inline function div_uc(i, j, k, grid, advection, U, c)
    1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, advective_tracer_flux_x, advection, U.u, c) +
                             δyᵃᶜᵃ(i, j, k, grid, advective_tracer_flux_y, advection, U.v, c) +
                             δzᵃᵃᶜ(i, j, k, grid, advective_tracer_flux_z, advection, U.w, c))
end


text:
@testset "warlpiri-inverses" begin

    #roundsinf test          ->          0000            0001            0010            0011            0100            0101            0110            0111            1000            1001            1010            1011            1100            1101            1110            1111
    const warlpiri_add_invs = [(WP0000→WP0000) (WP1111→WP0000) (WP1110→WP0000) (WP1101→WP0000) (WP1100→WP0000) (WP1011→WP0000) (WP1010→WP0000) (WP1001→WP0000) (WP1000→WP0000) (WP0111→WP0000) (WP0110→WP0000) (WP0101→WP0000) (WP0100→WP0000) (WP0011→WP0000) (WP0010→WP0000) (WP0001→WP0000)  # 0000 ->
                               (WP0000→WP1111) (WP1111→WP1111) (WP1110→WP1111) (WP1101→WP1111) (WP1100→WP1111) (WP1011→WP1111) (WP1010→WP1111) (WP1001→WP1111) (WP1000→WP1111) (WP0111→WP1111) (WP0110→WP1111) (WP0101→WP1111) (WP0100→WP1111) (WP0011→WP1111) (WP0010→WP1111) (WP0001→WP1111)  # 0001 ->
                               (WP0000→WP1110) (WP1111→WP1110) (WP1110→WP1110) (WP1101→WP1110) (WP1100→WP1110) (WP1011→WP1110) (WP1010→WP1110) (WP1001→WP1110) (WP1000→WP1110) (WP0111→WP1110) (WP0110→WP1110) (WP0101→WP1110) (WP0100→WP1110) (WP0011→WP1110) (WP0010→WP1110) (WP0001→WP1110)  # 0010 ->
                               (WP0000→WP1101) (WP1111→WP1101) (WP1110→WP1101) (WP1101→WP1101) (WP1100→WP1101) (WP1011→WP1101) (WP1010→WP1101) (WP1001→WP1101) (WP1000→WP1101) (WP0111→WP1101) (WP0110→WP1101) (WP0101→WP1101) (WP0100→WP1101) (WP0011→WP1101) (WP0010→WP1101) (WP0001→WP1101)  # 0011 ->
                               (WP0000→WP1100) (WP1111→WP1100) (WP1110→WP1100) (WP1101→WP1100) (WP1100→WP1100) (WP1011→WP1100) (WP1010→WP1100) (WP1001→WP1100) (WP1000→WP1100) (WP0111→WP1100) (WP0110→WP1100) (WP0101→WP1100) (WP0100→WP1100) (WP0011→WP1100) (WP0010→WP1100) (WP0001→WP1100)  # 0100 ->
                               (WP0000→WP1011) (WP1111→WP1011) (WP1110→WP1011) (WP1101→WP1011) (WP1100→WP1011) (WP1011→WP1011) (WP1010→WP1011) (WP1001→WP1011) (WP1000→WP1011) (WP0111→WP1011) (WP0110→WP1011) (WP0101→WP1011) (WP0100→WP1011) (WP0011→WP1011) (WP0010→WP1011) (WP0001→WP1011)  # 0101 ->
                               (WP0000→WP1010) (WP1111→WP1010) (WP1110→WP1010) (WP1101→WP1010) (WP1100→WP1010) (WP1011→WP1010) (WP1010→WP1010) (WP1001→WP1010) (WP1000→WP1010) (WP0111→WP1010) (WP0110→WP1010) (WP0101→WP1010) (WP0100→WP1010) (WP0011→WP1010) (WP0010→WP1010) (WP0001→WP1010)  # 0110 ->
                               (WP0000→WP1001) (WP1111→WP1001) (WP1110→WP1001) (WP1101→WP1001) (WP1100→WP1001) (WP1011→WP1001) (WP1010→WP1001) (WP1001→WP1001) (WP1000→WP1001) (WP0111→WP1001) (WP0110→WP1001) (WP0101→WP1001) (WP0100→WP1001) (WP0011→WP1001) (WP0010→WP1001) (WP0001→WP1001)  # 0111 ->
                               (WP0000→WP1000) (WP1111→WP1000) (WP1110→WP1000) (WP1101→WP1000) (WP1100→WP1000) (WP1011→WP1000) (WP1010→WP1000) (WP1001→WP1000) (WP1000→WP1000) (WP0111→WP1000) (WP0110→WP1000) (WP0101→WP1000) (WP0100→WP1000) (WP0011→WP1000) (WP0010→WP1000) (WP0001→WP1000)  # 1000 ->
                               (WP0000→WP0111) (WP1111→WP0111) (WP1110→WP0111) (WP1101→WP0111) (WP1100→WP0111) (WP1011→WP0111) (WP1010→WP0111) (WP1001→WP0111) (WP1000→WP0111) (WP0111→WP0111) (WP0110→WP0111) (WP0101→WP0111) (WP0100→WP0111) (WP0011→WP0111) (WP0010→WP0111) (WP0001→WP0111)  # 1001 ->
                               (WP0000→WP0110) (WP1111→WP0110) (WP1110→WP0110) (WP1101→WP0110) (WP1100→WP0110) (WP1011→WP0110) (WP1010→WP0110) (WP1001→WP0110) (WP1000→WP0110) (WP0111→WP0110) (WP0110→WP0110) (WP0101→WP0110) (WP0100→WP0110) (WP0011→WP0110) (WP0010→WP0110) (WP0001→WP0110)  # 1010 ->
                               (WP0000→WP0101) (WP1111→WP0101) (WP1110→WP0101) (WP1101→WP0101) (WP1100→WP0101) (WP1011→WP0101) (WP1010→WP0101) (WP1001→WP0101) (WP1000→WP0101) (WP0111→WP0101) (WP0110→WP0101) (WP0101→WP0101) (WP0100→WP0101) (WP0011→WP0101) (WP0010→WP0101) (WP0001→WP0101)  # 1011 ->
                               (WP0000→WP0100) (WP1111→WP0100) (WP1110→WP0100) (WP1101→WP0100) (WP1100→WP0100) (WP1011→WP0100) (WP1010→WP0100) (WP1001→WP0100) (WP1000→WP0100) (WP0111→WP0100) (WP0110→WP0100) (WP0101→WP0100) (WP0100→WP0100) (WP0011→WP0100) (WP0010→WP0100) (WP0001→WP0100)  # 1100 ->
                               (WP0000→WP0011) (WP1111→WP0011) (WP1110→WP0011) (WP1101→WP0011) (WP1100→WP0011) (WP1011→WP0011) (WP1010→WP0011) (WP1001→WP0011) (WP1000→WP0011) (WP0111→WP0011) (WP0110→WP0011) (WP0101→WP0011) (WP0100→WP0011) (WP0011→WP0011) (WP0010→WP0011) (WP0001→WP0011)  # 1101 ->
                               (WP0000→WP0010) (WP1111→WP0010) (WP1110→WP0010) (WP1101→WP0010) (WP1100→WP0010) (WP1011→WP0010) (WP1010→WP0010) (WP1001→WP0010) (WP1000→WP0010) (WP0111→WP0010) (WP0110→WP0010) (WP0101→WP0010) (WP0100→WP0010) (WP0011→WP0010) (WP0010→WP0010) (WP0001→WP0010)  # 1110 ->
                               (WP0000→WP0001) (WP1111→WP0001) (WP1110→WP0001) (WP1101→WP0001) (WP1100→WP0001) (WP1011→WP0001) (WP1010→WP0001) (WP1001→WP0001) (WP1000→WP0001) (WP0111→WP0001) (WP0110→WP0001) (WP0101→WP0001) (WP0100→WP0001) (WP0011→WP0001) (WP0010→WP0001) (WP0001→WP0001)  # 1111 ->
    ]

    add_inv_matrix = map(-, wmatrix)

    for (m, n) in zip(add_inv_matrix, warlpiri_add_invs)
        @test m == n
    end

    #roundsinf test          ->          0000            0001            0010            0011            0100            0101            0110            0111            1000            1001            1010            1011            1100            1101            1110            1111
    const warlpiri_mul_invs = [(WP1000→WP1000) (WP0111→WP1000) (WP0110→WP1000) (WP0101→WP1000) (WP0100→WP1000) (WP0011→WP1000) (WP0010→WP1000) (WP0001→WP1000) (WP0000→WP1000) (WP1111→WP1000) (WP1110→WP1000) (WP1101→WP1000) (WP1100→WP1000) (WP1011→WP1000) (WP1010→WP1000) (WP0000→WP1111)  # 0000 ->
                               (WP0001→WP0000) (WP0111→WP0111) (WP0110→WP0111) (WP0101→WP0111) (WP0100→WP0111) (WP0011→WP0111) (WP0010→WP0111) (WP0001→WP0111) (WP0000→WP0111) (WP1111→WP0111) (WP1110→WP0111) (WP1101→WP0111) (WP1100→WP0111) (WP1011→WP0111) (WP1010→WP0111) (WP1001→WP0111)  # 0001 ->
                               (WP1000→WP0110) (WP0111→WP0110) (WP0110→WP0110) (WP0101→WP0110) (WP0100→WP0110) (WP0011→WP0110) (WP0010→WP0110) (WP0001→WP0110) (WP0000→WP0110) (WP1111→WP0110) (WP1110→WP0110) (WP1101→WP0110) (WP1100→WP0110) (WP1011→WP0110) (WP1010→WP0110) (WP1001→WP0110)  # 0010 ->
                               (WP1000→WP0101) (WP0111→WP0101) (WP0110→WP0101) (WP0101→WP0101) (WP0100→WP0101) (WP0011→WP0101) (WP0010→WP0101) (WP0001→WP0101) (WP0000→WP0101) (WP1111→WP0101) (WP1110→WP0101) (WP1101→WP0101) (WP1100→WP0101) (WP1011→WP0101) (WP1010→WP0101) (WP1001→WP0101)  # 0011 ->
                               (WP1000→WP0100) (WP0111→WP0100) (WP0110→WP0100) (WP0101→WP0100) (WP0100→WP0100) (WP0011→WP0100) (WP0010→WP0100) (WP0001→WP0100) (WP0000→WP0100) (WP1111→WP0100) (WP1110→WP0100) (WP1101→WP0100) (WP1100→WP0100) (WP1011→WP0100) (WP1010→WP0100) (WP1001→WP0100)  # 0100 ->
                               (WP1000→WP0011) (WP0111→WP0011) (WP0110→WP0011) (WP0101→WP0011) (WP0100→WP0011) (WP0011→WP0011) (WP0010→WP0011) (WP0001→WP0011) (WP0000→WP0011) (WP1111→WP0011) (WP1110→WP0011) (WP1101→WP0011) (WP1100→WP0011) (WP1011→WP0011) (WP1010→WP0011) (WP1001→WP0011)  # 0101 ->
                               (WP1000→WP0010) (WP0111→WP0010) (WP0110→WP0010) (WP0101→WP0010) (WP0100→WP0010) (WP0011→WP0010) (WP0010→WP0010) (WP0001→WP0010) (WP0000→WP0010) (WP1111→WP0010) (WP1110→WP0010) (WP1101→WP0010) (WP1100→WP0010) (WP1011→WP0010) (WP1010→WP0010) (WP1001→WP0010)  # 0110 ->
                               (WP1000→WP0001) (WP0111→WP0001) (WP0110→WP0001) (WP0101→WP0001) (WP0100→WP0001) (WP0011→WP0001) (WP0010→WP0001) (WP0001→WP0001) (WP0000→WP0001) (WP1111→WP0001) (WP1110→WP0001) (WP1101→WP0001) (WP1100→WP0001) (WP1011→WP0001) (WP1010→WP0001) (WP1001→WP0001)  # 0111 ->
                               (WP1000→WP0000) (WP0111→WP0000) (WP0110→WP0000) (WP0101→WP0000) (WP0100→WP0000) (WP0011→WP0000) (WP0010→WP0000) (WP1000→WP0111) (WP0000→WP0000) (WP1111→WP0000) (WP1110→WP0000) (WP1101→WP0000) (WP1100→WP0000) (WP1011→WP0000) (WP1010→WP0000) (WP1001→WP0000)  # 1000 ->
                               (WP1000→WP1111) (WP0111→WP1111) (WP0110→WP1111) (WP0101→WP1111) (WP0100→WP1111) (WP0011→WP1111) (WP0010→WP1111) (WP0001→WP1111) (WP1001→WP1000) (WP1111→WP1111) (WP1110→WP1111) (WP1101→WP1111) (WP1100→WP1111) (WP1011→WP1111) (WP1010→WP1111) (WP1001→WP1111)  # 1001 ->
                               (WP1000→WP1110) (WP0111→WP1110) (WP0110→WP1110) (WP0101→WP1110) (WP0100→WP1110) (WP0011→WP1110) (WP0010→WP1110) (WP0001→WP1110) (WP0000→WP1110) (WP1111→WP1110) (WP1110→WP1110) (WP1101→WP1110) (WP1100→WP1110) (WP1011→WP1110) (WP1010→WP1110) (WP1001→WP1110)  # 1010 ->
                               (WP1000→WP1101) (WP0111→WP1101) (WP0110→WP1101) (WP0101→WP1101) (WP0100→WP1101) (WP0011→WP1101) (WP0010→WP1101) (WP0001→WP1101) (WP0000→WP1101) (WP1111→WP1101) (WP1110→WP1101) (WP1101→WP1101) (WP1100→WP1101) (WP1011→WP1101) (WP1010→WP1101) (WP1001→WP1101)  # 1011 ->
                               (WP1000→WP1100) (WP0111→WP1100) (WP0110→WP1100) (WP0101→WP1100) (WP0100→WP1100) (WP0011→WP1100) (WP0010→WP1100) (WP0001→WP1100) (WP0000→WP1100) (WP1111→WP1100) (WP1110→WP1100) (WP1101→WP1100) (WP1100→WP1100) (WP1011→WP1100) (WP1010→WP1100) (WP1001→WP1100)  # 1100 ->
                               (WP1000→WP1011) (WP0111→WP1011) (WP0110→WP1011) (WP0101→WP1011) (WP0100→WP1011) (WP0011→WP1011) (WP0010→WP1011) (WP0001→WP1011) (WP0000→WP1011) (WP1111→WP1011) (WP1110→WP1011) (WP1101→WP1011) (WP1100→WP1011) (WP1011→WP1011) (WP1010→WP1011) (WP1001→WP1011)  # 1101 ->
                               (WP1000→WP1010) (WP0111→WP1010) (WP0110→WP1010) (WP0101→WP1010) (WP0100→WP1010) (WP0011→WP1010) (WP0010→WP1010) (WP0001→WP1010) (WP0000→WP1010) (WP1111→WP1010) (WP1110→WP1010) (WP1101→WP1010) (WP1100→WP1010) (WP1011→WP1010) (WP1010→WP1010) (WP1001→WP1010)  # 1110 ->
                               (WP1000→WP1001) (WP0111→WP1001) (WP0110→WP1001) (WP0101→WP1001) (WP0100→WP1001) (WP0011→WP1001) (WP0010→WP1001) (WP0001→WP1001) (WP0000→WP1001) (WP1111→WP1001) (WP1110→WP1001) (WP1101→WP1001) (WP1100→WP1001) (WP1011→WP1001) (WP1010→WP1001) (WP1001→WP1001)  # 1111 ->=#
    ]

    mul_inv_matrix = map(/, wmatrix)

    for (m, n, w) in zip(mul_inv_matrix, warlpiri_mul_invs, wmatrix)
        @test (m, w) == (n, w)
    end
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
module PowerSystem

using Printf
using SparseArrays

import Base: show, get

const PQ_BUS_TYPE = 1
const PV_BUS_TYPE = 2
const REF_BUS_TYPE  = 3

"""
    AbstractPowerSystem

First layer of the package. Store the topology of a given
transmission network, including:

- the power injection at each bus ;
- the admittance matrix ;
- the default voltage at each bus.

Data are imported either from a matpower file, or a PSSE file.

"""
abstract type AbstractPowerSystem end

"""
    AbstractNetworkElement

Abstraction for all physical elements being parts of a `AbstractPowerSystem`.
Elements are divided in

- transmission lines (`Lines`)
- buses (`Buses`)
- generators (`Generators`)

"""
abstract type AbstractNetworkElement end

"""
    Buses <: AbstractNetworkElement

Buses of a transmission network.
"""
struct Buses <: AbstractNetworkElement end

"""
    Lines <: AbstractNetworkElement

Lines of a transmission network.
"""
struct Lines <: AbstractNetworkElement end

"""
    Generators <: AbstractElement

Generators in a transmission network
"""
struct Generators <: AbstractNetworkElement end

"""
    AbstractNetworkAttribute

Attribute of a `AbstractPowerSystem`.

"""
abstract type AbstractNetworkAttribute end

"""
    NumberOfBuses <: AbstractNetworkAttribute

Number of buses in a `AbstractPowerSystem`.
"""
struct NumberOfBuses <: AbstractNetworkAttribute end

"""
    NumberOfLines <: AbstractNetworkAttribute

Number of lines in a `AbstractPowerSystem`.
"""
struct NumberOfLines <: AbstractNetworkAttribute end

"""
    NumberOfGenerators <: AbstractNetworkAttribute

Number of generators in a `AbstractPowerSystem`.
"""
struct NumberOfGenerators <: AbstractNetworkAttribute end

"""
    NumberOfPVBuses <: AbstractNetworkAttribute

Number of PV buses in a `AbstractPowerSystem`.
"""
struct NumberOfPVBuses <: AbstractNetworkAttribute end

"""
    NumberOfPQBuses <: AbstractNetworkAttribute

Number of PQ buses in a `AbstractPowerSystem`.
"""
struct NumberOfPQBuses <: AbstractNetworkAttribute end

"""
    NumberOfSlackBuses <: AbstractNetworkAttribute

Number of slack buses in a `AbstractPowerSystem`.
"""
struct NumberOfSlackBuses <: AbstractNetworkAttribute end

"""
    BaseMVA <: AbstractNetworkAttribute

Base MVA of the network.
"""
struct BaseMVA <: AbstractNetworkAttribute end

"""
    BusAdmittanceMatrix <: AbstractNetworkAttribute

Bus admittance matrix associated with the topology of the network.
"""
struct BusAdmittanceMatrix <: AbstractNetworkAttribute end

abstract type AbstractIndexing <: AbstractNetworkAttribute end

"""
    PVIndexes <: AbstractIndexing

Indexes of the PV buses in a `AbstractPowerSystem`.
"""
struct PVIndexes <: AbstractIndexing end

"""
    PQIndexes <: AbstractIndexing

Indexes of the PQ buses in a `AbstractPowerSystem`.
"""
struct PQIndexes <: AbstractIndexing end

"""
    SlackIndexes <: AbstractIndexing

Indexes of the slack buses in a `AbstractPowerSystem`.
"""
struct SlackIndexes <: AbstractIndexing end

"""
    GeneratorIndexes <: AbstractIndexing

Indexes of the generators in a `AbstractPowerSystem`.
"""
struct GeneratorIndexes <: AbstractIndexing end

"""
    AbstractNetworkValues

Numerical values attached to the different attributes of the network.
"""
abstract type AbstractNetworkValues end

"""
    VoltageMagnitude <: AbstractNetworkValues

Magnitude `|v|` of the voltage `v = |v| exp(i θ)`.
"""
struct VoltageMagnitude <: AbstractNetworkValues end

"""
    VoltageAngle <: AbstractNetworkValues

Angle `θ` of the voltage `v = |v| exp(i θ)`.
"""
struct VoltageAngle <: AbstractNetworkValues end

"""
    ActivePower <: AbstractNetworkValues

Active power `P` of the complex power `S = P + iQ`.
"""
struct ActivePower <: AbstractNetworkValues end

"""
    ReactivePower <: AbstractNetworkValues

Reactive power `Q` of the complex power `S = P + iQ`.
"""
struct ReactivePower <: AbstractNetworkValues end

"""
    ActiveLoad <: AbstractNetworkValues

Active load `Pd` at buses.
"""
struct ActiveLoad <: AbstractNetworkValues end

"""
    ReactiveLoad <: AbstractNetworkValues

Reactive load `Qd` at buses.
"""
struct ReactiveLoad <: AbstractNetworkValues end

# Templating
"""
    get(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute)

Return value of attribute `attr` in the `AbstractPowerSystem` object `pf`.

    get(pf::AbstractPowerSystem, attr::AbstractIndexing)

Return indexing corresponding to a subset of the buses.

## Examples

```julia
npq = get(pf, NumberOfPQBuses())
npv = get(pf, NumberOfPVBuses())
index_pv = get(pf, PVIndexes()())
index_gen = get(pf, GeneratorIndexes()())

```
"""
function get end

"""
    bounds(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute, val::AbstractNetworkValues)

Return lower and upper bounds corresponding to the admissible values
of the `AbstractNetworkAttribute` `attr`.

## Examples

```julia
p_min, p_max = bounds(pf, Generator(), ActivePower())
v_min, v_max = bounds(pf, Buses(), VoltageMagnitude())

```
"""
function bounds end

struct Branches{T}
    Yff::Vector{T}
    Yft::Vector{T}
    Ytf::Vector{T}
    Ytt::Vector{T}
    from_buses::Vector{Int}
    to_buses::Vector{Int}
end

# Utils
function get_bus_id_to_indexes(bus)
    BUS_I = IndexSet.idx_bus()[1]
    nbus = size(bus, 1)

    bus_id_to_indexes = Dict{Int,Int}()
    for i in 1:nbus
        busn = bus[i, BUS_I]
        bus_id_to_indexes[busn] = i
    end
    return bus_id_to_indexes
end

function get_bus_generators(buses, gens, bus_id_to_indexes)
    bus_gen = Dict{Int, Array{Int}}()
    GEN_BUS = IndexSet.idx_gen()[1]
    for g in 1:size(gens, 1)
        bus_id = gens[g, GEN_BUS]
        bus_index = bus_id_to_indexes[bus_id]
        if haskey(bus_gen, bus_index)
            push!(bus_gen[bus_index], g)
        else
            bus_gen[bus_index] = Int[g]
        end
    end
    return bus_gen
end

function get_active_branches(lines, remove_lines)
    if isempty(remove_lines)
        return lines
    else
        lines = lines[1:end .!= remove_lines, :]
        return lines
    end
end

function import_dataset(datafile::String)
    if endswith(datafile, ".raw")
        data_raw = ParsePSSE.parse_raw(datafile)
        return ParsePSSE.raw_to_exapf(data_raw)
    elseif endswith(datafile, ".m")
        data_mat = ParseMAT.parse_mat(datafile)
        return ParseMAT.mat_to_exapf(data_mat)
    else
        error("Unsupported format in file $(datafile): supported extensions are " *
              "Matpower (.m) or PSSE (.raw)")
    end
end

include("indexes.jl")
using .IndexSet
include("parsers/parse_mat.jl")
using .ParseMAT
include("parsers/parse_psse.jl")
using .ParsePSSE

include("topology.jl")
include("power_network.jl")
include("matpower.jl")

end


text:
import LibGit2
import Pkg

# To manually generate the docs:
#
# 1. Install all python dependencies from the PYDEPS array below.
# 2. Run "julia make.jl" to generate the docs

# Set to true to disable some checks and cleanup
DEBUG = false

# Where to get files from and where to build them
SRCPATH = joinpath(@__DIR__, "src")
BUILDPATH = joinpath(@__DIR__, "build")
ROOTPATH = joinpath(@__DIR__, "..")
CONTINUOUS_INTEGRATION = get(ENV, "CI", nothing) == "true"

# Python and Julia dependencies needed for running the notebooks
PYDEPS = ["ase", "pymatgen"]
JLDEPS = [
    Pkg.PackageSpec(url="https://github.com/JuliaMolSim/DFTK.jl.git",
                    rev=LibGit2.head(ROOTPATH)),  # The current DFTK
]

# Setup julia dependencies for docs generation if not yet done
if !isfile(joinpath(@__DIR__, "Manifest.toml"))
    Pkg.activate(@__DIR__)
    Pkg.develop(Pkg.PackageSpec(path=ROOTPATH))
    Pkg.instantiate()
end

# Import packages for docs generation
using DFTK
using Documenter
using Literate

# Collect examples from the example index (src/index.md)
# The chosen examples are taken from the examples/ folder to be processed by Literate
EXAMPLES = [String(m[1])
            for m in match.(r"\"(examples/[^\"]+.md)\"",
                            readlines(joinpath(SRCPATH, "index.md")))
            if !isnothing(m)]

# Collect files to treat with Literate (i.e. the examples and the .jl files in the docs)
# The examples go to docs/literate_build/examples, the .jl files stay where they are
literate_files = [(src=joinpath(ROOTPATH, splitext(file)[1] * ".jl"),
                   dest=joinpath(SRCPATH, "examples"), example=true)
                  for file in EXAMPLES]
for (dir, directories, files) in walkdir(SRCPATH)
    for file in files
        if endswith(file, ".jl")
            push!(literate_files, (src=joinpath(dir, file), dest=dir, example=false))
        end
    end
end

# Function to insert badges to examples
function add_badges(str)
    badges = [
        "[![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/examples/@__NAME__.ipynb)",
        "[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/@__NAME__.ipynb)",
    ]

    # Find the Header and insert the badges right below
    splitted = split(str, "\n")
    idx = findfirst(startswith.(splitted, "# # "))
    idx === nothing && error("Example files must start with # #")
    insert!(splitted, idx + 1, "#md # " * badges[1])
    insert!(splitted, idx + 2, "#md # " * badges[2])
    insert!(splitted, idx + 3, "#md #")
    join(splitted, "\n")
end

# Run Literate on them all
for file in literate_files
    preprocess = file.example ? add_badges : identity
    Literate.markdown(file.src, file.dest; documenter=true, credit=false,
                      preprocess=preprocess)
    Literate.notebook(file.src, file.dest; credit=false,
                      execute=CONTINUOUS_INTEGRATION || DEBUG)
end

# Generate the docs in BUILDPATH
makedocs(
    modules = [DFTK],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = CONTINUOUS_INTEGRATION,
        canonical = "https://juliamolsim.github.io/DFTK.jl/stable/",
        assets = ["assets/favicon.ico"],
    ),
    sitename = "DFTK.jl",
    authors = "Michael F. Herbst, Antoine Levitt and contributors.",
    linkcheck = false,  # TODO
    linkcheck_ignore = [
        # Ignore links that point to GitHub's edit pages, as they redirect to the
        # login screen and cause a warning:
        r"https://github.com/([A-Za-z0-9_.-]+)/([A-Za-z0-9_.-]+)/edit(.*)",
    ],
    pages = [
        "Home" => "index.md",
        "Getting started" => Any[
            "guide/installation.md",
            "Tutorial" => "guide/tutorial.md",
            "guide/parallelization.md",
            "Density-functional theory" => "guide/density_functional_theory.md",
        ],
        "Examples" => EXAMPLES,
        "Advanced topics" => Any[
            "advanced/conventions.md",
            "advanced/data_structures.md",
            "advanced/useful_formulas.md",
            "advanced/symmetries.md",
        ],
        "api.md",
        "publications.md",
    ],
    strict = !DEBUG,
)

# Dump files for managing dependencies in binder
if CONTINUOUS_INTEGRATION
    cd(BUILDPATH) do
        open("environment.yml", "w") do io
            print(io,
                    """
                    name: dftk
                    channels:
                      - defaults
                      - conda-forge
                    dependencies:
                    """)
            for dep in PYDEPS
                println(io, "  - " * dep)
            end
        end

        # Install Julia dependencies into build
        Pkg.activate(".")
        for dep in JLDEPS
            Pkg.add(dep)
        end
    end
    Pkg.activate(@__DIR__)  # Back to Literate / Documenter environment
end

# Deploy docs to gh-pages branch
deploydocs(
    repo = "github.com/JuliaMolSim/DFTK.jl.git",
)

# Remove generated example files
if !DEBUG
    for file in literate_files
        base = splitext(basename(file.src))[1]
        for ext in [".ipynb", ".md"]
            rm(joinpath(file.dest, base * ext), force=true)
        end
    end
end

if !CONTINUOUS_INTEGRATION
    println("\nDocs generated, try $(joinpath(BUILDPATH, "index.html"))")
end


text:
function checkmatrix(m)
    size(m,1) != size(m,2) && error("Only square matrices of Numbers accepted, sorry!!!");
end


text:
const PyJuliaModuleValue_Type__ref = Ref(PyNULL)
PyJuliaModuleValue_Type() = begin
    ptr = PyJuliaModuleValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.ModuleValue",
                base = base,
                methods = [
                    (name = "seval", flags = Py_METH_VARARGS, meth = pyjlmodule_seval),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaModuleValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaModuleValue_New(x::Module) = PyJuliaValue_New(PyJuliaModuleValue_Type(), x)
PyJuliaValue_From(x::Module) = PyJuliaModuleValue_New(x)

pyjlmodule_seval(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsBetween("seval", args, 1, 2)) && return PyNULL
    if PyTuple_Size(args) == 1
        m = PyJuliaValue_GetValue(xo)::Module
        ism1(PyArg_GetArg(String, "seval", args, 0)) && return PyNULL
        s = takeresult(String)
    else
        ism1(PyArg_GetArg(Module, "seval", args, 0)) && return PyNULL
        m = takeresult(Module)
        ism1(PyArg_GetArg(String, "seval", args, 1)) && return PyNULL
        s = takeresult(String)
    end
    try
        PyObject_From(Base.eval(m, Meta.parse(s)))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end