text:
@testset "mechanism modification" begin
    @testset "flip_direction" begin
        Bounds = RigidBodyDynamics.Bounds
        jointname = "joint1"
        axis = [1.0, 0.0, 0.0]
        for JT in [Revolute, Prismatic]
            jointtype = JT(axis)
            posbounds = [Bounds(1., 2.)]
            velbounds = [Bounds(3., 4.)]
            effbounds = [Bounds(5., 6.)]
            joint = Joint(jointname, jointtype, position_bounds = posbounds, velocity_bounds = velbounds, effort_bounds = effbounds)
            flipped = RigidBodyDynamics.Graphs.flip_direction(joint)
            @test flipped.name == joint.name
            @test joint_type(flipped) == JT(.-axis)
            @test flipped.position_bounds == [Bounds(-2., -1.)]
            @test flipped.velocity_bounds == [Bounds(-4., -3.)]
            @test flipped.effort_bounds == [Bounds(-6., -5.)]
        end
    end

    @testset "attach!" begin
        Random.seed!(47)
        body0 = RigidBody{Float64}("root")
        mechanism = Mechanism(body0)
        joint1 = Joint("joint1", rand(Revolute{Float64}))
        joint1_pose = rand(Transform3D{Float64}, frame_before(joint1), default_frame(body0))
        body1 = RigidBody(rand(SpatialInertia{Float64}, frame_after(joint1)))
        joint2 = Joint("joint2", QuaternionFloating{Float64}())
        joint2_pose = rand(Transform3D{Float64}, frame_before(joint2), default_frame(body1))
        body2 = RigidBody(rand(SpatialInertia{Float64}, CartesianFrame3D("2")))

        # can't attach if predecessor is not among bodies of mechanism
        @test_throws AssertionError attach!(mechanism, body1, body2, joint2, joint_pose = joint2_pose)

        # attach body1
        attach!(mechanism, body0, body1, joint1, joint_pose = joint1_pose)
        @test length(bodies(mechanism)) == 2
        @test body1 ∈ bodies(mechanism)
        @test length(joints(mechanism)) == 1
        @test joint1 ∈ joints(mechanism)
        @test isapprox(fixed_transform(mechanism, frame_before(joint1), joint1_pose.to), joint1_pose)

        # can't use the same joint twice
        @test_throws AssertionError attach!(mechanism, body0, body1, joint1, joint_pose = joint1_pose)

        # attach body2
        attach!(mechanism, body1, body2, joint2, joint_pose = joint2_pose)
        @test length(bodies(mechanism)) == 3
        @test body2 ∈ bodies(mechanism)
        @test length(joints(mechanism)) == 2
        @test joint2 ∈ joints(mechanism)
        @test isapprox(fixed_transform(mechanism, frame_before(joint2), joint2_pose.to), joint2_pose)
    end

    @testset "attach! mechanism" begin
        Random.seed!(48)
        mechanism = randmech()
        nq = num_positions(mechanism)
        nv = num_velocities(mechanism)

        mechanism2 = rand_tree_mechanism(Float64, [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 5]; QuaternionSpherical{Float64}; Planar{Float64}; [Prismatic{Float64} for i = 1 : 5]]...)
        additional_frames = Dict{CartesianFrame3D, RigidBody{Float64}}()
        for body in bodies(mechanism2)
            for i = 1 : 5
                frame = CartesianFrame3D("frame_$i")
                tf = rand(Transform3D, frame, default_frame(body))
                add_frame!(body, tf)
                additional_frames[frame] = body
            end
        end

        parent_body = rand(collect(bodies(mechanism)))
        attach!(mechanism, parent_body, mechanism2)

        @test num_positions(mechanism) == nq + num_positions(mechanism2)
        @test num_velocities(mechanism) == nv + num_velocities(mechanism2)

        # make sure all of the frame definitions got copied over
        for frame in keys(additional_frames)
            body = additional_frames[frame]
            if body == root_body(mechanism2)
                body = parent_body
            end
            @test RigidBodyDynamics.is_fixed_to_body(body, frame)
        end

        state = MechanismState(mechanism) # issue 63
        rand!(state)
        M = mass_matrix(state)

        # independent acrobots in the same configuration
        # make sure mass matrix is block diagonal, and that blocks on diagonal are the same
        urdf = joinpath(@__DIR__, "urdf", "Acrobot.urdf")
        double_acrobot = parse_urdf(urdf)
        acrobot2 = parse_urdf(urdf)
        xsingle = MechanismState(acrobot2)
        rand!(xsingle)
        qsingle = configuration(xsingle)
        nq_single = length(qsingle)
        parent_body = root_body(double_acrobot)
        attach!(double_acrobot, parent_body, acrobot2)
        x = MechanismState(double_acrobot)
        set_configuration!(x, [qsingle; qsingle])
        H = mass_matrix(x)
        H11 = H[1 : nq_single, 1 : nq_single]
        H12 = H[1 : nq_single, nq_single + 1 : end]
        H21 = H[nq_single + 1 : end, 1 : nq_single]
        H22 = H[nq_single + 1 : end, nq_single + 1 : end]
        @test isapprox(H11, H22)
        @test isapprox(H12, zero(H12))
        @test isapprox(H21, zero(H21))
    end

    @testset "remove fixed joints" begin
        Random.seed!(49)
        joint_types = [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 10]; QuaternionSpherical{Float64}; Planar{Float64}; [Fixed{Float64} for i = 1 : 10]; [SinCosRevolute{Float64} for i = 1 : 5]]
        shuffle!(joint_types)
        mechanism = rand_tree_mechanism(Float64, joint_types...)
        state = MechanismState(mechanism)
        rand!(state)
        q = configuration(state)
        M = mass_matrix(state)
        nonfixedjoints = collect(filter(j -> !(joint_type(j) isa Fixed), tree_joints(mechanism)))

        remove_fixed_tree_joints!(mechanism)
        @test tree_joints(mechanism) == nonfixedjoints
        state_no_fixed_joints = MechanismState(mechanism)
        set_configuration!(state_no_fixed_joints, q)
        M_no_fixed_joints = mass_matrix(state_no_fixed_joints)
        @test isapprox(M_no_fixed_joints, M, atol = 1e-12)
    end

    @testset "replace joint" begin
        Random.seed!(50)
        joint_types = [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 10]; QuaternionSpherical{Float64}; Planar{Float64}; [Fixed{Float64} for i = 1 : 10]; [SinCosRevolute{Float64} for i = 1 : 5]]
        shuffle!(joint_types)
        mechanism = rand_tree_mechanism(Float64, joint_types...)

        for m in [mechanism; maximal_coordinates(mechanism)]
            for i = 1 : 10
                oldjoint = rand(joints(mechanism))
                newjoint = Joint("new", frame_before(oldjoint), frame_after(oldjoint), Fixed{Float64}())
                replace_joint!(mechanism, oldjoint, newjoint)
                @test newjoint ∈ joints(mechanism)
                @test oldjoint ∉ joints(mechanism)
            end
        end
    end

    @testset "submechanism" begin
        Random.seed!(51)
        for testnum = 1 : 100
            # joint_types = [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 10]; QuaternionSpherical{Float64}; Planar{Float64}; [Fixed{Float64} for i = 1 : 10]] # FIXME: use this
            joint_types = [Revolute{Float64} for i = 1 : 4]
            shuffle!(joint_types)
            mechanism = rand_tree_mechanism(Float64, joint_types...)
            state = MechanismState(mechanism)
            rand!(state)
            M = mass_matrix(state)

            submechanism_root = rand(collect(bodies(mechanism)))
            mechanism_part, bodymap, jointmap = submechanism(mechanism, submechanism_root)
            @test root_body(mechanism_part) == bodymap[submechanism_root]
            @test mechanism.gravitational_acceleration.v == mechanism_part.gravitational_acceleration.v

            substate = MechanismState(mechanism_part)
            for (oldjoint, newjoint) in jointmap
                set_configuration!(substate, newjoint, configuration(state, oldjoint))
                set_velocity!(substate, newjoint, velocity(state, oldjoint))
            end
            Msub = mass_matrix(substate)
            if num_velocities(mechanism_part) > 0
                indices = vcat([velocity_range(state, joint) for joint in tree_joints(mechanism) if joint ∈ keys(jointmap)]...)
                @test isapprox(M[indices, indices], Msub, atol = 1e-10)
            end
        end
    end

    @testset "reattach" begin
        Random.seed!(52)
        for testnum = 1 : 10
            # create random floating mechanism with flippable joints
            joint_types = [[Prismatic{Float64} for i = 1 : 10]; [Revolute{Float64} for i = 1 : 10]; [Fixed{Float64} for i = 1 : 10]; [SinCosRevolute{Float64} for i = 1 : 5]]
            shuffle!(joint_types)
            mechanism1 = rand_floating_tree_mechanism(Float64, joint_types...)

            # random state
            x1 = MechanismState(mechanism1)
            rand!(x1)

            # copy and set up mapping from bodies and joints of mechanism1 to those of mechanism2
            bodies1 = collect(bodies(mechanism1))
            joints1 = collect(joints(mechanism1))
            (bodies2, joints2, mechanism2) = deepcopy((bodies1, joints1, mechanism1))
            bodymap = Dict(zip(bodies1, bodies2))
            jointmap = Dict(zip(joints1, joints2))

            # find world, floating joint, floating body of mechanism1, and determine a new floating body
            world = root_body(mechanism1)
            floatingjoint = first(out_joints(world, mechanism1))
            floatingbody = successor(floatingjoint, mechanism1)
            newfloatingbody = rand(collect(non_root_bodies(mechanism1)))

            # reroot mechanism2
            newfloatingjoint = Joint("new_floating", QuaternionFloating{Float64}())
            joint_to_world = one(Transform3D, frame_before(newfloatingjoint), default_frame(world))
            body_to_joint = one(Transform3D, default_frame(newfloatingbody), frame_after(newfloatingjoint))
            attach!(mechanism2, bodymap[world], bodymap[newfloatingbody], newfloatingjoint, joint_pose = joint_to_world, successor_pose = body_to_joint)
            flipped_joint_map = Dict()
            remove_joint!(mechanism2, jointmap[floatingjoint]; flipped_joint_map = flipped_joint_map)

            # mimic the same state for the rerooted mechanism
            # copy non-floating joint configurations and velocities
            x2 = MechanismState(mechanism2)
            for (joint1, joint2) in jointmap
                if joint1 != floatingjoint
                    joint2_rerooted = get(flipped_joint_map, joint2, joint2)
                    set_configuration!(x2, joint2_rerooted, configuration(x1, joint1))
                    set_velocity!(x2, joint2_rerooted, velocity(x1, joint1))
                end
            end

            # set configuration and velocity of new floating joint
            newfloatingjoint_transform = inv(joint_to_world) * relative_transform(x1, body_to_joint.from, joint_to_world.to) * inv(body_to_joint)
            set_configuration!(x2, newfloatingjoint, newfloatingjoint_transform)

            newfloatingjoint_twist = transform(x1, relative_twist(x1, newfloatingbody, world), body_to_joint.from)
            newfloatingjoint_twist = transform(newfloatingjoint_twist, body_to_joint)
            newfloatingjoint_twist = Twist(body_to_joint.to, joint_to_world.from, newfloatingjoint_twist.frame, angular(newfloatingjoint_twist), linear(newfloatingjoint_twist))
            set_velocity!(velocity(x2, newfloatingjoint), newfloatingjoint, newfloatingjoint_twist)

            # do dynamics and compute spatial accelerations
            result1 = DynamicsResult(mechanism1)
            dynamics!(result1, x1)
            spatial_accelerations!(result1, x1)
            result2 = DynamicsResult(mechanism2)
            dynamics!(result2, x2)
            spatial_accelerations!(result2, x2)

            # make sure that joint accelerations for non-floating joints are the same
            for (joint1, joint2) in jointmap
                if joint1 != floatingjoint
                    joint2_rerooted = get(flipped_joint_map, joint2, joint2)
                    v̇1 = view(result1.v̇, velocity_range(x1, joint1))
                    v̇2 = view(result2.v̇, velocity_range(x2, joint2_rerooted))
                    @test isapprox(v̇1, v̇2)
                end
            end

            # make sure that body spatial accelerations are the same
            for (body1, body2) in bodymap
                accel1 = relative_acceleration(result1, body1, world)
                accel2 = relative_acceleration(result2, body2, bodymap[world])
                @test isapprox(angular(accel1), angular(accel2))
                @test isapprox(linear(accel1), linear(accel2))
            end
        end # for
    end # reattach

    @testset "maximal coordinates" begin
        Random.seed!(53)
        # create random tree mechanism and equivalent mechanism in maximal coordinates
        joint_types = [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 10]; QuaternionSpherical{Float64}; Planar{Float64}; [Fixed{Float64} for i = 1 : 10]; [SinCosRevolute{Float64} for i = 1 : 5]]
        tree_mechanism = rand_tree_mechanism(Float64, joint_types...)
        mc_mechanism, newfloatingjoints, bodymap, jointmap = maximal_coordinates(tree_mechanism)

        # randomize state of tree mechanism
        tree_state = MechanismState(tree_mechanism)
        rand!(tree_state)

        # put maximal coordinate system in state that is equivalent to tree mechanism state
        mc_state = MechanismState(mc_mechanism)
        for oldbody in non_root_bodies(tree_mechanism)
            newbody = bodymap[oldbody]
            joint = newfloatingjoints[newbody]

            tf = relative_transform(tree_state, frame_after(joint), frame_before(joint))
            set_configuration!(configuration(mc_state, joint), joint, tf)

            twist = transform(relative_twist(tree_state, frame_after(joint), frame_before(joint)), inv(tf))
            set_velocity!(mc_state, joint, twist)
        end
        setdirty!(mc_state)

        # do dynamics and compute spatial accelerations
        tree_dynamics_result = DynamicsResult(tree_mechanism);
        dynamics!(tree_dynamics_result, tree_state)
        spatial_accelerations!(tree_dynamics_result, tree_state)

        mc_dynamics_result = DynamicsResult(mc_mechanism);
        dynamics!(mc_dynamics_result, mc_state)
        spatial_accelerations!(mc_dynamics_result, mc_state)

        # compare spatial accelerations of bodies
        for (treebody, mcbody) in bodymap
            tree_accel = relative_acceleration(tree_dynamics_result, treebody, root_body(tree_mechanism))
            mc_accel = relative_acceleration(mc_dynamics_result, mcbody, root_body(mc_mechanism))
            @test isapprox(tree_accel, mc_accel; atol = 1e-10)
        end
    end # maximal coordinates

    @testset "generic scalar dynamics" begin # TODO: move to a better place
        Random.seed!(54)
        joint_types = [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 10]; QuaternionSpherical{Float64}; Planar{Float64}; [Fixed{Float64} for i = 1 : 10]; [SinCosRevolute{Float64} for i = 1 : 5]]
        mechanism = rand_tree_mechanism(Float64, joint_types...);
        mechanism, newfloatingjoints, bodymap, jointmap = maximal_coordinates(mechanism)

        state_float64 = MechanismState(mechanism)
        rand!(state_float64)
        NullDual = typeof(ForwardDiff.Dual(0., ()))
        state_dual = MechanismState{NullDual}(mechanism)
        configuration(state_dual) .= configuration(state_float64)
        velocity(state_dual) .= velocity(state_float64)

        dynamics_result_float64 = DynamicsResult(mechanism)
        dynamics!(dynamics_result_float64, state_float64)

        dynamics_result_dual = DynamicsResult{NullDual}(mechanism)
        dynamics!(dynamics_result_dual, state_dual)

        @test isapprox(dynamics_result_float64.v̇, dynamics_result_dual.v̇; atol = 1e-3)
        @test isapprox(dynamics_result_float64.λ, dynamics_result_dual.λ; atol = 1e-3)
    end

    @testset "modcount" begin
        Random.seed!(55)
        joint_types = [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1 : 10]; QuaternionSpherical{Float64}; Planar{Float64}; [Fixed{Float64} for i = 1 : 10]; [SinCosRevolute{Float64} for i = 1 : 5]]
        mechanism = rand_tree_mechanism(Float64, joint_types...);
        state = MechanismState(mechanism)
        attach!(mechanism, root_body(mechanism), RigidBody(rand(SpatialInertia{Float64}, CartesianFrame3D())), Joint("bla", QuaternionFloating{Float64}()))
        @test_throws RigidBodyDynamics.ModificationCountMismatch mass_matrix(state)
        state2 = MechanismState(mechanism)
        mass_matrix(state2) # make sure this doesn't throw
        @test_throws RigidBodyDynamics.ModificationCountMismatch mass_matrix(state) # make sure this still throws
        try
            mass_matrix(state)
        catch e
            showerror(devnull, e)
        end
    end
end # mechanism modification

ext: jl


text:
using ADCME
using AdFem

include("../chip_unstructured_solver.jl")
include("../chip_unstructured_geometry.jl")

k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

function k_exact(x, y)
    k_mold + 1000 * k_chip_ref * (x-0.49)^2 / (1 + x^2)
end

function k_nn(x, y, θ)
    out =  fc(x, [20,20,20,1], θ)^2  .+ k_mold
    squeeze(out)
end

xy = mesh.nodes
xy2 = zeros(mesh.nedge, 2)
for i = 1:mesh.nedge
    xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
end
xy = [xy;xy2]

x, y = xy[chip_fem_idx, 1], xy[chip_fem_idx, 2]

n = 901
θ = Variable(randn(n))
k_chip = k_nn(x, y, θ)
k_chip_exact = @. k_exact(x, y)

loss =  mean((k_chip .- k_chip_exact)^2)
loss = loss * 1e10

# ---------------------------------------------------
_loss = Float64[]
cb = (vs, iter, loss)->begin
    global _loss
    push!(_loss, loss)
    printstyled("[#iter $iter] loss=$loss\n", color=:green)
end

ADCME.options.training.training = placeholder(true)
l = placeholder(rand(193,))
x = placeholder(rand(5875, 2))

# train the neural network
opt = AdamOptimizer().minimize(loss)
sess = Session(); init(sess)
for i = 1:10
    _, loss_ = run(sess, [opt, loss], feed_dict=Dict(l=>k_chip_exact, x=>xy))
    @info i, loss_
end

# sess = Session(); init(sess)
BFGS!(sess, loss, vars = [k_chip], callback = cb)
matwrite("k_diagnose_theta.mat", Dict("theta"=>run(sess, θ)))
matwrite("k_diagnose_loss.mat", Dict("loss"=>_loss))


text:
"""
    mlm_backest_sum_perms(data; nPerms, permFun, isXIntercept, isZIntercept,
                          weights, targetType, isXSum, isZSum, isMainEff)

Obtains permutation p-values for MLM t-statistics, including back-estimated
left-out sum contrast interactions.

# Arguments

- data = RawData object
- nPerms = number of permutations. Defaults to `1000`.

# Keyword arguments

- permFun = function used to permute `Y`. Defaults to `shuffle_rows`
  (shuffles rows of `Y`).
- isXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- isZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- weights = 1d array of floats to use as column weights for `Y`, or `nothing`.
  If the former, must be the same length as the number of columns of `Y`.
  Defaults to `nothing`.
- targetType = string indicating the target type toward which to shrink the
  error variance, or `nothing`. If the former, acceptable inputs are "A", "B",
  "C", and "D". Defaults to `nothing`.
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries
- isMainEff = boolean flag indicating whether or not to include p-values for
  the main effects
- isXSum = boolean flag indicating whether an left-out `X` sum contrast
  interaction needs to be back-estimated
- isZSum = boolean flag indicating whether an left-out `Z` sum contrast
  interaction needs to be back-estimated

# Value

Tuple
- tStats: 2d array of floats; t-statistics
- pvals: 2d array of floats; permutation p-values

# Some notes

Permutations are computed in parallel when possible.

This implementation assumes that the `X` predictor matrix
(when `isXSum` is set to `true`) and the `Z` predictor matrix
(when `isZSum` is set to `true`) consist of the sum contrasts for a single
categorical variable (and no other variables).

"""
function mlm_backest_sum_perms(data::RawData, nPerms::Int64=1000;
                               permFun::Function=shuffle_rows,
                               isXIntercept::Bool=true,
                               isZIntercept::Bool=true,
                               weights=nothing, targetType=nothing,
                               isXSum::Bool=true, isZSum::Bool=true,
                               isMainEff::Bool=false)

    # Wrapper function that performs MLM, back-estimates left-out sum sum
    # contrast interactions, and gets t-statistics
    function mlm_backest_sum_t_stat(data::RawData)
        return t_stat(mlm_backest_sum(data; isXIntercept=isXIntercept,
                                      isZIntercept=isZIntercept,
                                      weights=weights, targetType=targetType,
                                      isXSum=isXSum, isZSum=isZSum), isMainEff)
    end

    # Get MLM t-statistcs and permutation p-values
    return perm_pvals(mlm_backest_sum_t_stat, data, nPerms; permFun=permFun)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -31.730360167469108 <= q <= 25.84895899492546)

                     add_NL_constraint(m, :(softplus(0.0048766473277672695 + 0.6370542772950833*softplus(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + 0.7379502189935647*softplus(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + -0.017781498731094914*softplus(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + -0.30964053180650675*softplus(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + -0.19693223852539354*softplus(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + softplus(0.6750746912592471 + -0.6689814871888631*softplus(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + 0.08174772592079638*softplus(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + 0.43040165086652493*softplus(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + 0.6360324545158482*softplus(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + 0.2590621611090005*softplus(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + softplus(-0.5603004507305789 + -0.6080760485894823*softplus(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + -0.17827205367530619*softplus(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + -0.3024929990799108*softplus(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + 0.2876553916986979*softplus(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + 0.1394466763243627*softplus(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + softplus(-0.3104079613918116 + -0.8967757209559282*softplus(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + 0.9554659562201504*softplus(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + 0.4284452595185697*softplus(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + -0.8546417933790886*softplus(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + -0.9550224035691546*softplus(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + softplus(0.11524778075403264 + -0.08842997013200282*softplus(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + -0.2446780712327863*softplus(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + 0.5457935761049009*softplus(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + -0.1673283141057702*softplus(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + 0.6343389483483679*softplus(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("algorithm.jl");

numberOfParticle = 10;
numberOfOperands = 3;
minRang = 100;
maxRang = 200;
targetValue = 100;
maxIteration = 2000
maxVelocity = 10

swarm = generateSwarm(numberOfParticle,numberOfOperands,minRang,maxRang);
PsoAlgorithm(swarm,targetValue,maxIteration);

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using Documenter
using Literate
using vSmartMOM
using vSmartMOM.Absorption
using vSmartMOM.Scattering
using vSmartMOM.vSmartMOM
using Plots

ENV["GKSwstype"] = "nul"

function build()

    tutorials = ["Tutorial_Absorption.jl", "Tutorial_Scattering.jl"] # ,
    tutorials_paths = [joinpath(@__DIR__, "src", "pages", "tutorials", tutorial) for tutorial in tutorials]

    for tutorial in tutorials_paths
        Literate.markdown(tutorial, joinpath(@__DIR__, "src", "pages", "tutorials"))
    end

    tutorials_md = [joinpath("pages", "tutorials", tutorial[1:end-3]) * ".md" for tutorial in tutorials]

    pages = Any[
        "Getting Started"       => "index.md",
        "vSmartMOM"             => Any[
                                    "Overview" => "pages/vSmartMOM/Overview.md",
                                    "Example" => "pages/vSmartMOM/Example.md",
                                    "User-Defined RT Parameters" => "pages/vSmartMOM/InputParametersGuide.md",
                                    "Methods & Types" => "pages/vSmartMOM/Types.md",
                                    "References" => "pages/vSmartMOM/References.md"
                                    ],
        "Absorption"            => [
                                    "Overview" => "pages/Absorption/Overview.md",
                                    "Example" => "pages/Absorption/Example.md",
                                    "Methods & Types" => "pages/Absorption/Types.md",
                                    "References" => "pages/Absorption/References.md"
                                    ],
        "Scattering"            => [
                                    "Overview" => "pages/Scattering/Overview.md",
                                    "Example" => "pages/Scattering/Example.md",
                                    "Methods & Types" => "pages/Scattering/Types.md",
                                    "References" => "pages/Scattering/References.md"
                                    ],
        "Tutorials"             => tutorials_md
    ]

    mathengine = MathJax(Dict(
        :TeX => Dict(
            :equationNumbers => Dict(:autoNumber => "AMS"),
            :Macros => Dict(),
        ),
    ))

    # The format will make other pages in parallel with the index page
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = mathengine,
        collapselevel = 1,
        assets = ["assets/favicon.ico"]
    )

    # This way it shows warnings of functions that have not been documented
    makedocs(
        sitename = "vSmartMOM",
        format = format,
        clean = false,
        modules = [vSmartMOM],
        pages = pages,
    )

end

build()

deploydocs(
    repo = "github.com/RemoteSensingTools/vSmartMOM.jl.git",
    target = "build",
    push_preview = true,
)


text:
# Device type and auxiliary functions

import Base: convert, @deprecate_binding

export
    devcount,
    CuDevice, name, totalmem, attribute, capability,
    list_devices

typealias CuDevice_t Cint

"Return the number of available CUDA devices"
function devcount()
    count_ref = Ref{Cint}()
    @apicall(:cuDeviceGetCount, (Ptr{Cint},), count_ref)
    return count_ref[]
end

immutable CuDevice
    ordinal::Cint
    handle::CuDevice_t

    function CuDevice(i::Integer)
        ordinal = convert(Cint, i)
        handle_ref = Ref{CuDevice_t}()
        @apicall(:cuDeviceGet, (Ptr{CuDevice_t}, Cint), handle_ref, ordinal)
        new(ordinal, handle_ref[])
    end
end

convert(::Type{CuDevice_t}, dev::CuDevice) = dev.handle

"Get the name of a CUDA device"
function name(dev::CuDevice)
    const buflen = 256
    buf = Array(Cchar, buflen)
    @apicall(:cuDeviceGetName, (Ptr{Cchar}, Cint, CuDevice_t),
                               buf, buflen, dev)
    buf[end] = 0
    return unsafe_string(pointer(buf))
end

"Get the amount of GPU memory (in bytes) of a CUDA device"
function totalmem(dev::CuDevice)
    mem_ref = Ref{Csize_t}()
    @apicall(:cuDeviceTotalMem, (Ptr{Csize_t}, CuDevice_t), mem_ref, dev)
    return mem_ref[]
end

@enum(CUdevice_attribute, MAX_THREADS_PER_BLOCK = Cint(1),
                          MAX_BLOCK_DIM_X,
                          MAX_BLOCK_DIM_Y,
                          MAX_BLOCK_DIM_Z,
                          MAX_GRID_DIM_X,
                          MAX_GRID_DIM_Y,
                          MAX_GRID_DIM_Z,
                          MAX_SHARED_MEMORY_PER_BLOCK,
                          TOTAL_CONSTANT_MEMORY,
                          WARP_SIZE,
                          MAX_PITCH,
                          MAX_REGISTERS_PER_BLOCK,
                          CLOCK_RATE,
                          TEXTURE_ALIGNMENT,
                          GPU_OVERLAP,
                          MULTIPROCESSOR_COUNT,
                          KERNEL_EXEC_TIMEOUT,
                          INTEGRATED,
                          CAN_MAP_HOST_MEMORY,
                          COMPUTE_MODE,
                          MAXIMUM_TEXTURE1D_WIDTH,
                          MAXIMUM_TEXTURE2D_WIDTH,
                          MAXIMUM_TEXTURE2D_HEIGHT,
                          MAXIMUM_TEXTURE3D_WIDTH,
                          MAXIMUM_TEXTURE3D_HEIGHT,
                          MAXIMUM_TEXTURE3D_DEPTH,
                          MAXIMUM_TEXTURE2D_LAYERED_WIDTH,
                          MAXIMUM_TEXTURE2D_LAYERED_HEIGHT,
                          MAXIMUM_TEXTURE2D_LAYERED_LAYERS,
                          SURFACE_ALIGNMENT,
                          CONCURRENT_KERNELS,
                          ECC_ENABLED,
                          PCI_BUS_ID,
                          PCI_DEVICE_ID,
                          TCC_DRIVER,
                          MEMORY_CLOCK_RATE,
                          GLOBAL_MEMORY_BUS_WIDTH,
                          L2_CACHE_SIZE,
                          MAX_THREADS_PER_MULTIPROCESSOR,
                          ASYNC_ENGINE_COUNT,
                          UNIFIED_ADDRESSING,
                          MAXIMUM_TEXTURE1D_LAYERED_WIDTH,
                          MAXIMUM_TEXTURE1D_LAYERED_LAYERS,
                          CAN_TEX2D_GATHER,
                          MAXIMUM_TEXTURE2D_GATHER_WIDTH,
                          MAXIMUM_TEXTURE2D_GATHER_HEIGHT,
                          MAXIMUM_TEXTURE3D_WIDTH_ALTERNATE,
                          MAXIMUM_TEXTURE3D_HEIGHT_ALTERNATE,
                          MAXIMUM_TEXTURE3D_DEPTH_ALTERNATE,
                          PCI_DOMAIN_ID,
                          TEXTURE_PITCH_ALIGNMENT,
                          MAXIMUM_TEXTURECUBEMAP_WIDTH,
                          MAXIMUM_TEXTURECUBEMAP_LAYERED_WIDTH,
                          MAXIMUM_TEXTURECUBEMAP_LAYERED_LAYERS,
                          MAXIMUM_SURFACE1D_WIDTH,
                          MAXIMUM_SURFACE2D_WIDTH,
                          MAXIMUM_SURFACE2D_HEIGHT,
                          MAXIMUM_SURFACE3D_WIDTH,
                          MAXIMUM_SURFACE3D_HEIGHT,
                          MAXIMUM_SURFACE3D_DEPTH,
                          MAXIMUM_SURFACE1D_LAYERED_WIDTH,
                          MAXIMUM_SURFACE1D_LAYERED_LAYERS,
                          MAXIMUM_SURFACE2D_LAYERED_WIDTH,
                          MAXIMUM_SURFACE2D_LAYERED_HEIGHT,
                          MAXIMUM_SURFACE2D_LAYERED_LAYERS,
                          MAXIMUM_SURFACECUBEMAP_WIDTH,
                          MAXIMUM_SURFACECUBEMAP_LAYERED_WIDTH,
                          MAXIMUM_SURFACECUBEMAP_LAYERED_LAYERS,
                          MAXIMUM_TEXTURE1D_LINEAR_WIDTH,
                          MAXIMUM_TEXTURE2D_LINEAR_WIDTH,
                          MAXIMUM_TEXTURE2D_LINEAR_HEIGHT,
                          MAXIMUM_TEXTURE2D_LINEAR_PITCH,
                          MAXIMUM_TEXTURE2D_MIPMAPPED_WIDTH,
                          MAXIMUM_TEXTURE2D_MIPMAPPED_HEIGHT,
                          COMPUTE_CAPABILITY_MAJOR,
                          COMPUTE_CAPABILITY_MINOR,
                          MAXIMUM_TEXTURE1D_MIPMAPPED_WIDTH,
                          STREAM_PRIORITIES_SUPPORTED,
                          GLOBAL_L1_CACHE_SUPPORTED,
                          LOCAL_L1_CACHE_SUPPORTED,
                          MAX_SHARED_MEMORY_PER_MULTIPROCESSOR,
                          MAX_REGISTERS_PER_MULTIPROCESSOR,
                          MANAGED_MEMORY,
                          MULTI_GPU_BOARD,
                          MULTI_GPU_BOARD_GROUP_ID)
@assert Cint(MULTI_GPU_BOARD_GROUP_ID) == 85
@deprecate_binding SHARED_MEMORY_PER_BLOCK MAX_SHARED_MEMORY_PER_BLOCK
@deprecate_binding REGISTERS_PER_BLOCK MAX_REGISTERS_PER_BLOCK
Base.deprecate(:GPU_OVERLAP)
@deprecate_binding MAXIMUM_TEXTURE2D_ARRAY_WIDTH MAXIMUM_TEXTURE2D_LAYERED_WIDTH
@deprecate_binding MAXIMUM_TEXTURE2D_ARRAY_HEIGHT MAXIMUM_TEXTURE2D_LAYERED_HEIGHT
@deprecate_binding MAXIMUM_TEXTURE2D_ARRAY_NUMSLICES MAXIMUM_TEXTURE2D_LAYERED_LAYERS
Base.deprecate(:CAN_TEX2D_GATHER)

function attribute(dev::CuDevice, attrcode::CUdevice_attribute)
    value_ref = Ref{Cint}()
    @apicall(:cuDeviceGetAttribute, (Ptr{Cint}, Cint, CuDevice_t),
                                    value_ref, attrcode, dev)
    return value_ref[]
end

"Return the compute capabilities of a CUDA device"
function capability(dev::CuDevice)
    major_ref = Ref{Cint}()
    minor_ref = Ref{Cint}()
    @apicall(:cuDeviceComputeCapability, (Ptr{Cint}, Ptr{Cint}, CuDevice_t),
                                         major_ref, minor_ref, dev)
    return VersionNumber(major_ref[], minor_ref[])
end

"List all CUDA devices with their capabilities and attributes"
function list_devices()
    cnt = devcount()
    if cnt == 0
        println("No CUDA-capable device found.")
        return
    end

    for i = 0:cnt-1
        dev = CuDevice(i)
        nam = name(dev)
        tmem = round(Integer, totalmem(dev) / (1024^2))
        cap = capability(dev)

        println("device[$i]: $(nam), capability $(cap.major).$(cap.minor), total mem = $tmem MB")
    end
end


text:
struct ShrikeIndex{T <: AbstractFloat}
    data::Array{T, 2}
    npoints::Int
    ndims::Int
    depth::Int
    ntrees::Int
    random_vectors::Array{T, 2}
    splits::Array{T, 2}
    indexes::Array{Array{Int,1}, 2}
end

"""
    ShrikeIndex(data::Array{T, 2}, depth::Int, ntrees::Int) where T -> ensemble

Constructor for ensemble of sparse random projection trees with voting. Returns `ShrikeIndex` type.
(An ensemble of multiple random projection trees.)

** Type Fields**
* `data::Array{T, 2}`: Contains all data points
* `npoints::Int`: Number of data points
* `ndims::Int`: Dimensionality of the data
* `depth::Int`: maximum depth of the tree. (Depth of 0 means only a root, depth of 1 means root has two children)
* `ntrees::Int`: Number of trees to make
* `random_vectors::AbstractArray`: The random projections used to make the tree
* `splits::Array{T, 2}`: The split values for each node in each tree stored in a 2D array
* `indexes::Array{Array{Int,1}, 2}`: 2D array of datapoint indexes at each leaf node in each tree.
Note that RP forest does not store indexes at non-leaf nodes.

Follows the implementation outlined in:

>**Fast Nearest Neighbor Search through Sparse Random Projections and Voting.**
>Ville Hyvönen, Teemu Pitkänen, Sotirios Tasoulis, Elias Jääsaari, Risto Tuomainen,
>Liang Wang, Jukka Ilmari Corander, Teemu Roos. Proceedings of the 2016 IEEE
>Conference on Big Data (2016)

with some modifications.
"""
function ShrikeIndex(data::AbstractArray{T, 2}, depth::Int, ntrees::Int) where T
    ndims, npoints = size(data)
    maxdepth = floor(Int, log(2.0, npoints) - log(2.0, 1))
    if depth > maxdepth
        depth = maxdepth
    end
    sparse_vecs::Bool = ndims > 500
    random_vectors = random_projections(T, depth*ntrees, ndims, sparse = sparse_vecs)
    nsplits = 2^depth -1 # Number of splits in the tree (Same as number of non-leaf nodes)
    splits = zeros(T, ntrees, nsplits) # Array for storing hyperplane split values
    projections = transpose(data) * random_vectors  # projection[i, j] = the projection of point i onto the jth random vector
    nleafs = 2^depth
    leaf_node_data_idxs = Array{Array{Int, 1}, 2}(undef, nleafs, ntrees) # Array for storing data indexes at the leaf

    for t in 1:ntrees
        prev_level = [collect(1:npoints)]
        next_level = Array{Array{Int, 1}, 2}(undef, 2, 1)
        next_split_idx = 1
        for d in 1:depth
            proj_row = (t-1)*depth + d # Current row of projection values
            for i in 1:length(prev_level)
                idxs = prev_level[i]
                # Compute median value of the projection of all current datapoints
                p = reshape(projections[idxs, proj_row], :)
                med::T = median(p) # I think this requires a sort, so this process can be sped up slightly
                splits[t, next_split_idx] = med
                next_split_idx += 1
                # Separate data points across hyperplane using the median
                mask = p .< med
                # Send points with the projection less than the median to the left child
                next_level[2*i - 1] = idxs[mask]
                # Send points with the projection greater than the median to the right child
                next_level[2*i] = idxs[.!mask]

            end
            prev_level = next_level
            next_level = Array{Array{Int, 1}, 2}(undef, 2^(d+1), 1)
        end
        leaf_node_data_idxs[:, t] = prev_level
    end

    return ShrikeIndex(data, npoints, ndims, depth, ntrees, random_vectors, splits, leaf_node_data_idxs)
end

ShrikeIndex(data::AbstractArray{T, 2}; depth::Int=4, ntrees::Int=5) where T = ShrikeIndex(data, depth, ntrees)

"""
    ShrikeIndex(data::AbstractArray{T, 2}, max)k; depth::Union{Int, Float64}=Inf, ntrees::Int=5) -> shi

Keyword argument version of the constructor that includes intended number of nearest neighbors.

If the default `depth` is used, the constructor sets the tree depth as deep as
possible given `max_k`. This way, the accuracy/memory tradeoff is determined directly by `ntrees` and the
desired `vote_cutoff` (`vote_cutoff` is a parameter passed to `ann` or `knngraph`).

If an argument is passed for `depth`, constructor attempts to use the supplied `depth`, but guarentees
that the depth of the tree is shallow enough to ensure that each leaf has at least k points.
(Without this check, the index may return less than k neighbors when queried.)

**Parameters**

1. `data`: A (dxn) array. Each column is a datapoint with dimension `d`.
2. `max_k`: The maximum number of neighbors that will be queried. If intend
to use the `ShrikeIndex` to approximate at most 10 nearest neigbors of a point,
set `max_k = 10`. This argument is used to infer the deepest tree depth possible
so as to maximize speed,

**Keyword Arguments**

3. `ntrees`: The number of trees in the index. More trees means more accuracy,
more memory and less speed. Use this to tune the speed/accuracy tradeoff.
4. `depth`: The number of splits in the tree. Depth of 0 means only a root,
depth of 1 means root has two children, etc..
"""
function ShrikeIndex(data::AbstractArray{T, 2}, maxk::Int; depth::Union{Int, Float64}=Inf, ntrees::Int=10) where T
    m, n = size(data)
    maxdepth = floor(Int, log(2, n) - log(2, maxk))
    if depth > maxdepth
        safedepth = maxdepth
    else
        safedepth = depth
    end
    return ShrikeIndex(data, safedepth, ntrees)
end

function Base.show(io::IO, shi::ShrikeIndex)
    descr = "ShrikeIndex: \n    $(shi.ntrees) trees \n    $(shi.npoints) datapoints \n    Depth $(shi.depth)"
    print(io, descr)
end

function random_projections(T::Type, nvecs::Int, ndims::Int; sparse::Bool=true)
    # Sparse arrays are not implemented because sparse x dense matmul
    # was slower than dense x dense matmul when benchmarked.
    # Additionally, I'm not sure how to use sparse arrays and still
    # be type stable.
    # a::Float64 = 1. / √ ndims
    # if sparse
    #     R = sprand(ndims, nvecs, a)
    # else
        R = randn(T, ndims, nvecs)
        R ./= sum(R.^2, dims=2) .^ 0.5
    # end
    return R::Array{T, 2}
end


text:
import LazySets.require

function __init__()
    @require JuMP = "4076af6c-e467-56ae-b986-b466b2749572" include(raw"initialization/jump.jl")
    @require MathematicalSystems = "d14a8603-c872-5ed3-9ece-53e0e82e39da" include(raw"initialization/mathematical_systems.jl")
    # This packages are already required in LazySets (i think)
    # @require Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" ""
    # @require CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
    # @require Optim = "429524aa-4258-5aef-a3af-852621145aeb"
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
using Test
using NastranResults

import NastranResults.readmeta!,
       NastranResults.readcomponents!,
       NastranResults.PCHLL,
       NastranResults.collect,
       NastranResults.PunchModalData,
       NastranResults.PMData,
       NastranResults.PunchFrequencyData,
       NastranResults.PFData,
       NastranResults.Component

# need to implement this soon
function unit_test_readop4_nofile()
    @test true
end

"""
    test readop4 constants
"""
function unit_test_op4_constants()
    @test NastranResults.COL == 65536
    @test NastranResults.iCOL == 1/NastranResults.COL
    @test NastranResults.MAT_ALL == "ALL"
    @test NastranResults.MAT_NAMES == "NAMES"
end

function unit_test_readop4_names()
    @test true
end

function unit_test_readop4_nomatrix()
    @test true
end

function unit_test_readop4_MAA()
    @test true
end

function unit_test_readop4_KAA()
    @test true
end

function unit_test_readop4_ALL()
    @test true
end

function unit_test_readop4_correct()
    @test true
end

function unit_test_readpch_nofile()
    @test true
end

"""
    test readpunch constants
"""
function unit_test_readpch_constants()
    @test NastranResults.TITLE == "TITLE"
    @test NastranResults.SUBTITLE == "SUBTITLE"
    @test NastranResults.LABEL == "LABEL"
    @test NastranResults.SUBCASE_ID == "SUBCASE ID"
    @test NastranResults.EIGENVALUE == "EIGENVALUE"
    @test NastranResults.MODE == "MODE"
    @test NastranResults.POINT == "POINT ID"
    @test NastranResults.INCLUDE == "INCLUDE"
    @test NastranResults.SINGLEP == 7
    @test NastranResults.DOUBLEP == 15
    @test NastranResults.PCHL == 18
    @test NastranResults.PCHLL == 72
    @test NastranResults.DOLLAR == '\$'
    @test NastranResults.CONT == "-CONT-"
    @test NastranResults.EQUAL == '='
    @test NastranResults.DOF == 6
end

"""
    test function readmeta! for modal type
"""
function unit_test_readpch_readmeta_modal(T = Float64)
    lines = []
    push!(lines, "\$TITLE   = EIGENMODES                                                          1")
    push!(lines, "\$SUBTITLE= SUBTITLE                                                            2")
    push!(lines, "\$LABEL   =EIGENMODES ANALYSIS                                                  3")
    push!(lines, "\$EIGENVECTOR                                                                   4")
    push!(lines, "\$REAL OUTPUT                                                                   5")
    push!(lines, "\$SUBCASE ID =           1                                                      6")
    push!(lines, "\$EIGENVALUE =  1.0000000E+00  MODE =     1                                     7")

    data = PunchModalData{T}()
    for line in lines
        line = line[1:PCHLL]
        readmeta!(data, line)
    end

    @test data.title == "EIGENMODES"
    @test data.subtitle == "SUBTITLE"
    @test data.label == "EIGENMODES ANALYSIS"
    @test data.datatype == "EIGENVECTOR"
    @test data.outputtype == "REAL OUTPUT"
    @test data.subcaseid == 1
    @test data.eigenvalue == 1.0
    @test data.mode == 1
end

"""
    test function reameta! for frequency type
"""
function unit_test_readpch_readmeta_frequency(T = Float64)
    lines = []
    push!(lines, "\$TITLE   = TITLE                                                               1")
    push!(lines, "\$SUBTITLE= 1 X                                                                 2")
    push!(lines, "\$LABEL   = 1 X                                                                 3")
    push!(lines, "\$DISPLACEMENTS                                                                 4")
    push!(lines, "\$REAL-IMAGINARY OUTPUT                                                         5")
    push!(lines, "\$SUBCASE ID =           2                                                      6")
    push!(lines, "\$POINT ID =           1  IDENTIFIED BY FREQUENCY                               7")

    data = PunchFrequencyData{T}()
    for line in lines
        line = line[1:PCHLL]
        readmeta!(data, line)
    end

    @test data.title == "TITLE"
    @test data.subtitle == "1 X"
    @test data.label == "1 X"
    @test data.datatype == "DISPLACEMENTS"
    @test data.outputtype == "REAL-IMAGINARY OUTPUT"
    @test data.subcaseid == 2
    @test data.point == 1
end

"""
    test function readcomponents! for modal type
"""
function unit_test_readpch_readcomponents_modal(T = Float64)
    l1 = "  00000001       G      1.000000E+00      2.000000E+00      3.000000E+00       1"
    l2 = "-CONT-                 -4.000000E+00      5.000000E+00      6.000000E+00       2"
    lines = [l1, l2]

    data = PunchModalData{T}()
    node = Vector{Int}(undef, 0)
    next = iterate(lines)
    (line, state) = next
    readcomponents!(data, node, line, lines, state)

    @test node[1] == 1
    @test data.eigenvector[1].x ==  1.0
    @test data.eigenvector[1].y ==  2.0
    @test data.eigenvector[1].z ==  3.0
    @test data.eigenvector[1].rx == -4.0
    @test data.eigenvector[1].ry ==  5.0
    @test data.eigenvector[1].rz ==  6.0
end

"""
    test function readcomponents! for frequency type
"""
function unit_test_readpch_readcomponents_frequency(T = Float64)
    l1 = "    0.000000E+00 G     -1.000000E+00      2.000000E+00      3.000000E+00       1"
    l2 = "-CONT-                 -4.000000E+00      5.000000E+00      6.000000E+00       2"
    l3 = "-CONT-                 -7.000000E+00      8.000000E+00      9.000000E+00       3"
    l4 = "-CONT-                 -1.000000E+01      1.100000E+01      1.200000E+01       4"
    lines = [l1, l2, l3, l4]

    data = PunchFrequencyData{T}()
    frequency = Vector{T}(undef, 0)
    next = iterate(lines)
    (line, state) = next
    readcomponents!(data, frequency, line, lines, state)

    @test frequency[1] == 0.0
    @test data.response[1].x == -1.0 - 7.0im
    @test data.response[1].y == 2.0 + 8.0im
    @test data.response[1].z == 3.0 + 9.0im
    @test data.response[1].rx == -4.0 - 10.0im
    @test data.response[1].ry == 5.0 + 11.0im
    @test data.response[1].rz == 6.0 + 12.0im
end

"""
    test function readpch! for frequency type
"""
function integration_test_readpch_modal(T = Float64)
    data = readpch("$(@__DIR__)/assets/modal.pch", PunchModalData{T})

    @test typeof(first(data.node)) == Int
    @test typeof(first(data.modal[1].eigenvector).x) == T
    @test data.node == [1]
    @test data.modal[1].title == "EIGENMODES"
    @test data.modal[1].subtitle == "SUBTITLE"
    @test data.modal[1].label == "EIGENMODES ANALYSIS"
    @test data.modal[1].datatype == "EIGENVECTOR"
    @test data.modal[1].outputtype == "REAL OUTPUT"
    @test data.modal[1].subcaseid == 1
    @test data.modal[1].eigenvalue == 1.0
    @test data.modal[1].mode ==  1
    @test data.modal[1].eigenvector[1].x == 1.0
    @test data.modal[1].eigenvector[1].y == 2.0
    @test data.modal[1].eigenvector[1].z == 3.0
    @test data.modal[1].eigenvector[1].rx == -4.0
    @test data.modal[1].eigenvector[1].ry == 5.0
    @test data.modal[1].eigenvector[1].rz == 6.0
end

"""
    test function readpch! for frequency type
"""
function integration_test_readpch_frequency(T = Float64)
    data = readpch("$(@__DIR__)/assets/frequency.pch", PunchFrequencyData{T})

    @test typeof(first(data.frequency)) == T
    @test typeof(first(data.response[1].response).x) == Complex{T}
    @test data.frequency == [0.0]
    @test data.response[1].title == "TITLE"
    @test data.response[1].subtitle == "1 X"
    @test data.response[1].label == "1 X"
    @test data.response[1].datatype == "DISPLACEMENTS"
    @test data.response[1].outputtype == "REAL-IMAGINARY OUTPUT"
    @test data.response[1].subcaseid == 2
    @test data.response[1].point ==  1
    @test data.response[1].response[1].x == -1.0 - 7.0im
    @test data.response[1].response[1].y == 2.0 + 8.0im
    @test data.response[1].response[1].z == 3.0 + 9.0im
    @test data.response[1].response[1].rx == -4.0 - 10.0im
    @test data.response[1].response[1].ry == 5.0 + 11.0im
    @test data.response[1].response[1].rz == 6.0 + 12.0im
end

"""
    test function collect for modal type
"""
function unit_test_collect_component(T = Float64)
    d = []
    n = 3
    testvector = randn(6, n)
    for i = 1:n
        push!(d, collect(Component{T}(testvector[:, i]...)))
    end

    for i in 1:n
        @test d[i] == testvector[:, i]
    end
end

"""
    test function collect for modal type
"""
function unit_test_collect_modal(T = Float64)
    data = Vector{PunchModalData{T}}(undef, 0)
    node = Vector{Int}(undef, 0)
    n = 3
    testvector = randn(6, n)
    for i = 1:n
        mdata = PunchModalData{T}()
        push!(mdata.eigenvector, Component{T}(testvector[:, i]...))
        push!(node, i)
        push!(data, mdata)
    end
    testdata = PMData{T}(data, node)
    d = collect(testdata)

    for i in 1:n
        @test d[:, i] == testvector[:, i]
    end
end

"""
    test function collect for frequency type
"""
function unit_test_collect_frequency(T = Float64)
    data = Vector{PunchFrequencyData{T}}(undef, 0)
    node = Vector{Int}(undef, 0)
    n = 3
    testvector = randn(6, n) + randn(6, n)*im
    for i = 1:n
        mdata = PunchFrequencyData{T}()
        push!(mdata.response, Component{Complex{T}}(testvector[:, i]...))
        push!(node, i)
        push!(data, mdata)
    end
    testdata = PFData{T}(data, node)
    d = collect(testdata)

    for i in 1:n
        @test d[:, i] == testvector[:, i]
    end
end

# these tests look more like integration tests but ok for now
@testset "Test OP4  " begin
    unit_test_readop4_nofile()
    unit_test_op4_constants()

    # need small test file
    #unit_test_readop4_nonames()
    #unit_test_readop4_nomatrix()
    #unit_test_readop4_MAA()
    #unit_test_readop4_KAA()
    #unit_test_readop4_ALL()
    #unit_test_readop4_correct()
end

@testset "Test PUNCH" begin
    unit_test_readpch_nofile()
    unit_test_readpch_constants()
    unit_test_readpch_readmeta_modal()
    unit_test_readpch_readmeta_frequency()
    unit_test_readpch_readcomponents_modal()
    unit_test_readpch_readcomponents_frequency()

    unit_test_collect_component()
    unit_test_collect_modal()
    unit_test_collect_frequency()

    integration_test_readpch_modal()
    integration_test_readpch_frequency()
end

#@testset "Test HDF5 " begin
#    @test true
#end

max_stars_repo_name: mgcth/SDNastran

max_issues_repo_name: mgcth/SDNastran

max_forks_repo_name: mgcth/SDNastran


text:
include("A2C_utils.jl")
using Flux, Statistics, Random, Distributions, StatsFuns, Zygote

T = 100
Nhidden = 50
Nstates = 3
Nstate_rep = 3
Naction = 2
Nin = Naction+1+1+Nstates #2 actions, 1 rew, 1 time, 3 states
Nout = Naction+1 #action 1, action 2, value function
βv = 0.05f0
βe = 0.05f0
γ = 0.9f0

### specify environment function per task ###
function env(y, t, state, ps)
    ### input the agent output; output action, reward, and new input
    πt = exp.(y[1:2]) #probability of actions
    d = Binomial(1, Float64(πt[1])) #probability of action 1
    a = (2 - rand(d)) # 1 -> 1, 0 -> 2

    ahot = zeros(2)
    if Bool(state[1]) #in initial state; move to new state
        prew = 0 #no reward
        p1 = [0.8; 0.2][a] #probability of going to state 2/3 is 0.8/0.2 for action 1 and vice versa for action 2
        s_new = [0; 0; 0]
        s = 3-rand(Binomial(1, p1)) #0->3, 1->2
        s_new[s] = 1
        #println(s, s_new)
        t = 0
        if rand() < 0.025 ps = [ps[2]; ps[1]] end #flip reward
        ahot[a] = 1.
    else
        if Bool(state[2]) prew = ps[1] else prew = ps[2] end #action-independent reward probability
        s_new = [1; 0; 0] #back to initial state
        t = 1
    end

    rew = rand(Binomial(1, Float64(prew))) #draw reward probabilistically
    x = [ahot; rew; t; s_new] #input is action, reward, timestep, state
    return Int(a), Float32(rew), Float32.(x), Float32.(ps), Float32.(s_new)
end
Zygote.@nograd env #don't take gradients of this

### task-specific initialization ###
function initialize(ps, state)
    if maximum(ps) <= 0
        if rand() < 0.5
            ps = [0.9; 0.1]
        else
            ps = [0.1; 0.9]
        end
    end
    if maximum(state) <= 0
        state = [1; 0; 0]
    end
    return Float32.(ps), Float32.(state)
end
Zygote.@nograd initialize #don't take gradients of this

### task specific evaluation/progress function ####
function model_eval()
    #evaluation function; compute mean reward (random is 0.5; oracle is 0.74)
    means = []
    all_as = []
    for i = 1:50
        L, ys, rews, as, ps = run_episode(m, hidden = false)
        means = [means; 2*mean(rews)] #compute mean reward for second stage
        all_as = [all_as; mean(as)]
    end
    return mean(means), mean(all_as)
end

max_stars_repo_name: KrisJensen/metaRL

max_issues_repo_name: KrisJensen/metaRL

max_forks_repo_name: KrisJensen/metaRL


text:
# Conjugate

export Conjugate

"""
**Convex conjugate**

    Conjugate(f)

Returns the convex conjugate (also known as Fenchel conjugate, or Fenchel-Legendre transform) of function `f`, that is
```math
f^*(x) = \\sup_y \\{ \\langle y, x \\rangle - f(y) \\}.
```
"""
struct Conjugate{T <: ProximableFunction} <: ProximableFunction
  f::T
  function Conjugate{T}(f::T) where {T<: ProximableFunction}
    if is_convex(f) == false
      error("`f` must be convex")
    end
    new(f)
  end
end

is_prox_accurate(f::Conjugate) = is_prox_accurate(f.f)
is_convex(f::Conjugate) = true
is_cone(f::Conjugate) = is_cone(f.f) && is_convex(f.f)
is_smooth(f::Conjugate) = is_strongly_convex(f.f)
is_strongly_convex(f::Conjugate) = is_smooth(f.f)
is_quadratic(f::Conjugate) = is_strongly_convex(f.f) && is_generalized_quadratic(f.f)
is_generalized_quadratic(f::Conjugate) = is_quadratic(f.f)

fun_dom(f::Conjugate) = fun_dom(f.f)

Conjugate(f::T) where {T <: ProximableFunction} = Conjugate{T}(f)

# only prox! is provided here, call method would require being able to compute
# an element of the subdifferential of the conjugate

function prox!(y::AbstractArray{R}, g::Conjugate, x::AbstractArray{R}, gamma::R=one(R)) where R <: Real
  # Moreau identity
  v = prox!(y, g.f, x/gamma, one(R)/gamma)
  if is_set(g)
    v = zero(R)
  else
    v = dot(x,y) - gamma*dot(y,y) - v
  end
  for k in eachindex(y)
    y[k] = x[k] - gamma*y[k]
  end
  return v
end

# complex case, need to cast inner products to real

function prox!(y::AbstractArray{Complex{T}}, g::Conjugate, x::AbstractArray{Complex{T}}, gamma::R=one(R)) where {R <: Real, T <: RealOrComplex{R}}
  v = prox!(y, g.f, x/gamma, one(R)/gamma)
  if is_set(g)
    v = zero(R)
  else
    v = real(dot(x,y)) - gamma*real(dot(y,y)) - v
  end
  for k in eachindex(y)
    y[k] = x[k] - gamma*y[k]
  end
  return v
end

# naive implementation

function prox_naive(g::Conjugate, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  y, v = prox_naive(g.f, x/gamma, 1.0/gamma)
  return x - gamma*y, real(dot(x,y)) - gamma*real(dot(y,y)) - v
end


text:
using Revise
using ADCME
using AdFem
using PyPlot
using SparseArrays
using MAT
using Random
Random.seed!(1234)

close("all")

# grid setup
m = 20
n = 20
h = 1/n
bc = bcedge("all", m, n, h)
bd = bcnode("all", m, n, h)
bd_2d = [bd; (m+1)*(n+1) .+ bd]

# time step setup
t = 0;
t_final = 0.01;
NT = 20;
dt = t_final/NT;

# physical constants
ρ = 1.0
μ = Variable(10.0)
ν = μ / ρ

# pre-compute constant matrices
mass_mat = constant(compute_fem_mass_matrix1(m, n, h))
mass_mat_2d = constant(compute_fem_mass_matrix(m, n, h))
mass_mat_2d_bdry, _ = fem_impose_Dirichlet_boundary_condition(mass_mat_2d, bd, m, n, h) # 2D boundary condition
laplace_mat = constant(compute_fem_laplace_matrix1(m, n, h))
laplace_mat_2d = constant(compute_fem_laplace_matrix(m, n, h))
interact_mat = constant(compute_interaction_matrix(m, n, h))
tpfa_mat_bdry, _ = compute_fvm_tpfa_matrix(ones(m*n), bc, zeros(size(bc,1)), m, n, h)
tpfa_mat_bdry = constant(tpfa_mat_bdry)

# exact solutions
function u_exact(x,y,t)
    cos(2*pi*x) * sin(2*pi*y) * exp(-8*pi*pi*ν*t)
    # cos(x) * sin(y) * exp(-2*ν*t)
end

function v_exact(x,y,t)
    -sin(2*pi*x) * cos(2*pi*y) * exp(-8*pi*pi*ν*t)
    # -sin(x) * cos(y) * exp(-2*ν*t)
end

function p_exact(x, y, t, ρ)
    -ρ/4 * (cos(4*pi*x) + cos(4*pi*y)) * exp(-16*pi*pi*ν*t)
    # -ρ/4 * (cos(2*x1) + cos(2*x2)) * exp(-4*ν*t)
end

function step1(U, p0, Source = missing)
    Source = coalesce(Source, zeros(2*(m+1)*(n+1)))
    u0 = U[1:(m+1)*(n+1)]
    v0 = U[(m+1)*(n+1)+1:end]
    u0_gauss = fem_to_gauss_points(u0, m, n, h)
    v0_gauss = fem_to_gauss_points(v0, m, n, h)

    gradu = eval_grad_on_gauss_pts1(u0, m, n, h) # du/dx = gradu[:,1], du/dy = gradu[:,2]
    gradv = eval_grad_on_gauss_pts1(v0, m, n, h) # dv/dx = gradv[:,1], dv/dy = gradv[:,2]

    M1 = mass_mat
    M2 = compute_fem_mass_matrix1(gradu[:,1], m, n, h)
    M3 = compute_fem_advection_matrix1(u0_gauss, v0_gauss, m, n, h)
    M4 = laplace_mat
    A11 = 1/dt * M1 + M2 + M3 + ν * M4

    A12 = compute_fem_mass_matrix1(gradu[:,2], m, n, h)

    A21 = compute_fem_mass_matrix1(gradv[:,1], m, n, h)

    M2 = compute_fem_mass_matrix1(gradv[:,2], m, n, h)
    A22 = 1/dt * M1 + M2 + M3 + ν * M4    # M1, M3, M4 are same as A11

    A = [A11 A12
        A21 A22]

    grad_p = compute_interaction_term(p0, m, n, h) # weak form of [dp/dx; dp/dy] on fem points

    s1 = u0_gauss .* gradu[:,1] + v0_gauss .* gradu[:,2]
    s2 = u0_gauss .* gradv[:,1] + v0_gauss .* gradv[:,2]
    b3 = compute_fem_source_term(s1, s2, m, n, h)

    F = Source + 1/ρ * grad_p - ν * laplace_mat_2d * [u0;v0] - b3

    A, _ = fem_impose_Dirichlet_boundary_condition(A, bd, m, n, h)
    F = scatter_update(F, bd_2d, zeros(length(bd_2d)))

    sol = A \ F
    return sol
end

function step2(U_int)
    rhs = ρ / dt * interact_mat * U_int
    sol = tpfa_mat_bdry \ rhs
    return sol
end

function step3(U_int, dp)
    grad_dp = - compute_interaction_term(dp, m, n, h)
    rhs = mass_mat_2d * U_int - dt/ρ * grad_dp
    rhs = scatter_update(rhs, bd_2d, zeros(length(bd_2d)))
    sol = mass_mat_2d_bdry \ rhs
    return sol
end

# input: U length 2(m+1)(n+1)
# input: p length mn
function solve_ns_one_step(U, p)
    dU = step1(U, p)
    U_int = U + dU
    dp = step2(U_int)
    p_new = p + dp
    U_new = step3(U_int, dp)
    return U_new, p_new
end

function condition(i, velo_arr, p_arr)
    i <= NT + 1
end

function body(i, velo_arr, p_arr)
    velo = read(velo_arr, i-1)
    pres = read(p_arr, i-1)
    # op = tf.print("i=",i)
    # i = bind(i, op)
    velo_new, pres_new = solve_ns_one_step(velo, pres)
    velo_arr = write(velo_arr, i, velo_new)
    p_arr = write(p_arr, i, pres_new)
    return i+1, velo_arr, p_arr
end

# fem nodes
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = @.  u_exact(x, y, 0.0)
u0=stack(u0)
v0 = @.  v_exact(x, y, 0.0)
v0=stack(v0)
velo_arr = TensorArray(NT+1)
velo_arr = write(velo_arr, 1, [u0; v0])

# fvm nodes
xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @.  p_exact(x, y, 0.0, ρ)
p0=stack(p0)
p_arr = TensorArray(NT+1)
p_arr = write(p_arr, 1, p0)

i = constant(2, dtype=Int32)

_, velo, p = while_loop(condition, body, [i, velo_arr, p_arr])
velo = set_shape(stack(velo), (NT+1, 2*(m+1)*(n+1)))
p = set_shape(stack(p), (NT+1, m*n))

velo_data = matread("ns_implicit_param_data.mat")["V"]
sample_size=20
idx = rand(1:(m+1)*(n+1), sample_size)
idx = [idx; idx .+ (m+1)*(n+1)] # observe velocity in y direction at the same location
observed_data = velo_data[:,idx]

noise = true
noise_level = 0.1
if noise
    noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
    observed_data = observed_data .* noise_ratio
end

loss = mean((velo[:,idx] - observed_data)^2)

max_iter = 100
sess = Session(); init(sess)
loss_ = BFGS!(sess, loss, max_iter)
figure(); semilogy(loss_); title("loss vs. iteration")
savefig(string("loss_ns_viscosity",string(noise_level),".png"))


text:
using MAGMA

using Test, LinearAlgebra

@testset "random singular precision matrices" begin

    matrixToTest = rand(Cfloat, 2, 2)

    right_answer = svd(matrixToTest).S
    S = right_answer

    jobu = 'A'
    jobvt = 'A'

    ldu=2
    ldvt=2
    lwork=400
    success=magmaInit()

    result = gesvd!(jobu,jobvt,matrixToTest)

    s = result[2]

    magmaFinalize()

    diff = S .- s
    error_value = norm(diff)

    @test error_value < 1e-7

    if error_value >= 1e-7
        println("Unfortunately, the test failed.")
        println("Here is some maybe useful information:")
        println("the element_type is ", typeof(matrixToTest))
        println("the right answer = ", S)
        println("However, MAGMA got the answer = ", s)
        println("The info returned by MAGMA is: ", info[1])
    end
end


text:
using ArgParse
using Compat
import Compat: UTF8String, ASCIIString, view, readstring, foreach

include("CLI.jl")
include("commands.jl")

function main()
    parsed_args = parse_args(build_arg_table())
    command = parsed_args["%COMMAND%"]

    if command == "help"
        foreach(x -> println(x), help())
    elseif command == "align_many"
        align_many(parsed_args[command]["algorithm"],
        parsed_args[command]["gap_penalty"],
        parsed_args[command]["extension_penalty"],
        parsed_args[command]["sequences"],
        parsed_args[command]["score_mat"],
        parsed_args[command]["output_file"],
        parsed_args[command]["write_alignments"]
        )
    elseif command == "compare_pr"
        compare_pr(parsed_args[command]["algorithm"],
        parsed_args[command]["pos_seqs"],
        parsed_args[command]["neg_seqs"],
        parsed_args[command]["score_mat"],
        parsed_args[command]["gp_min"],
        parsed_args[command]["gp_max"],
        parsed_args[command]["output_file"],
        parsed_args[command]["tpr"])
    elseif command == "roc"
        bench_roc(parsed_args[command]["algorithm"],
        parsed_args[command]["pos_seqs"],
        parsed_args[command]["neg_seqs"],
        parsed_args[command]["score_mat"],
        parsed_args[command]["gap_penalty"],
        parsed_args[command]["extension_penalty"],
        parsed_args[command]["output_file"])
    else
        align(parsed_args[command]["algorithm"],
        parsed_args[command]["gap_penalty"],
        parsed_args[command]["extension_penalty"],
        parsed_args[command]["sequence_1"],
        parsed_args[command]["sequence_2"],
        parsed_args[command]["score_mat"],
        parsed_args[command]["output_file"]
        )
    end
end

# fire up simulation if run using command line
if !isinteractive()
    main()
end


text:
module MyUtils

export subtype_tree, project

function subtype_tree_rekursive(t::Type, depth::Int = 0)
    println(string("     "^depth, "|--> ", t))
    subs = subtypes(t)
    for t in subs
        subtype_tree_rekursive(t, depth+1)
    end
end

"""
    subtype_tree(t::Type)

Show for a given type `t` the supertype in the first line,
this type in the second line and all derived types below.

# Examples

```jldoctest
julia> subtype_tree(AbstractFloat)
Real
|--> AbstractFloat
     |--> BigFloat
     |--> Float16
     |--> Float32
     |--> Float64
```
"""
function subtype_tree(t::Type)
    println(supertype(t))
    subtype_tree_rekursive(t, 0)
end

"""
    project(x::Real, x1::Real, y1::Real, x2::Real, y2::Real)

Calculate for a `x` from the range `[x1,y1]` its value in the range `[x2,y2]`.

This includes `x`s outside of the input range and reverse.
See the Examples section below.

# Examples

```jldoctest
julia> project(5.0,4.0,6.0,3.0,8.0)
5.5
julia> project(2.0,0.0,1.0,0.0,1.0)
2.0
julia> project(0.3,1.0,0.0,0.0,1.0)
0.7
```
"""
@inline function project(x::Real, x1::Real, y1::Real, x2::Real, y2::Real)
    return (x-x1)/(y1-x1)*(y2-x2)+x2
end

end


text:
include("useful_functions.jl")

function solve_003()
    N = 600851475143
    while true
        k = get_smallest_divisor(N)
        N::Int64 = N/k
        !isprime(N) || break
    end
    return N
end


text:
using ArgParse
using Cambrian
using CartesianGeneticProgramming
using IICGP
using LinearAlgebra
import Cambrian.mutate

"""
    function generate_io(n::Int64=10)::Tuple{Array{Float64,2},Array{Float64,1}}

Function generating random inputs/outputs with a simple function mapping.
The function mapping is
    y = 0.5 * (x[1] + x[2] * x[3])
where `x[i]` is in `[0.0, 10.0]`.
Can be used to generate a simple dataset.
"""
function generate_io(n::Int64=1000)::Tuple{Array{Float64,2},Array{Float64,1}}
    inps = 10 * rand(Float64, (n, 3))
    return inps, 0.5 * (inps[:, 1] .+ inps[:, 2] .* inps[:, 3])
end

"""
    function fitness(ind::CGPInd, inps::Array{Float64,2}, outs::Array{Float64,1})::Array{Float64,1}

Fitness function for float-CGP test.
Fitness is calculated based on the L2-error prediction from the given
inputs/outputs dataset.
"""
function fitness(ind::CGPInd, inps::Array{Float64,2},
                 outs::Array{Float64,1})::Array{Float64,1}
    preds = zeros(Float64, (size(inps)[1]))
    for i in eachindex(inps[:,1])
        preds[i] = CartesianGeneticProgramming.process(ind, inps[i,:])[1]
    end
    return [-LinearAlgebra.norm(preds .- outs, 1) / length(preds)]
end

"""
    function fitness(ind::CGPInd, n::Int64=100)::Array{Float64,1}

Fitness function for float-CGP test.
Fitness is calculated based on the L2-error prediction from the generated
inputs/outputs dataset.
"""
function fitness(ind::CGPInd, n::Int64=1000)::Array{Float64,1}
    inps, outs = generate_io(n)
    fitness(ind, inps, outs)
end

# Read configuration file
s = ArgParseSettings()
@add_arg_table! s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/flcgp.yaml"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end
args = parse_args(ARGS, s)
n_in = 3  # Three floating numbers as input
n_out = 1  # Single output
cfg = CartesianGeneticProgramming.get_config(args["cfg"]; n_in=n_in, n_out=n_out)

#=
# Test I/O without evolution
inp, out = generate_io(1)
bar = CGPInd(cfg)
CartesianGeneticProgramming.set_inputs(bar, inp)
CartesianGeneticProgramming.process(bar)
=#

#=
# Test perfect ind
function generate_perfect_ind(cfg::NamedTuple)
    chromosome = [
        0.3, 0.1, 0.0, 0.0, 0.0, # xs
        0.5, 0.6, 0.0, 0.0, 0.0, # ys
        0.4, 0.1, 0.1, 0.1, 0.1, # fs
        0.6  # output
    ]
    return CGPInd(cfg, chromosome)
end
foo = generate_perfect_ind(cfg)
=#

if length(args["ind"]) > 0
    ind = CGPInd(cfg, read(args["ind"], String))
    ftn = fitness(ind, inps, outs)
    println("Fitness: ", ftn)
else
    # Define mutate and fit functions
    mutate(i::CGPInd) = goldman_mutate(cfg, i)
    fit(i::CGPInd) = fitness(i, 1000)
    # Create an evolution framework
    e = CGPEvolution(cfg, fit)
    # Run evolution
    run!(e)
end


text:
const TableauVSPARK = AbstractTableauSPARK{:vspark}
const ParametersVSPARK = AbstractParametersVSPARK{:vspark}

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for Variational systems.

This integrator solves the following system of equations for the internal stages,
```math
\begin{align}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 ,
\end{align}
```
with definitions
```math
\begin{align}
P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{align}
```
and update rule
```math
\begin{align}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) .
\end{align}
```
"""
struct IntegratorVSPARK{DT, TT, tabType,
                                ET <: IDAE{DT,TT},
                                PT <: ParametersVSPARK{DT,TT},
                                ST <: NonlinearSolver{DT},
                                IT <: InitialGuessPODE{DT,TT}, D, S, R} <: AbstractIntegratorVSPARK{DT, TT}
    equation::ET
    tableau::AbstractTableauSPARK{tabType,TT}

    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheSPARK{DT,TT,D,S,R}
end

function IntegratorVSPARK(equation::IDAE{DT,TT},
                         tableau::AbstractTableauSPARK{ST,TT}, Δt::TT) where {DT,TT,ST}
    D = equation.d
    S = tableau.s
    R = tableau.r
    P = tableau.ρ

    @assert tableau.ρ == tableau.r-1

    N = 3*D*S + 3*D*R

    if isdefined(tableau, :d) && length(tableau.d) > 0
        N += D
    end

    # create params
    params = ParametersVSPARK{DT,D,S,R,P}(equation.f, equation.ϑ, equation.u, equation.g, equation.ϕ, Δt, tableau)

    # create solver
    solver = create_nonlinear_solver(DT, N, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheSPARK{DT,TT,D,S,R}()

    # create integrator
    IntegratorVSPARK{DT, TT, ST, typeof(equation), typeof(params), typeof(solver), typeof(iguess), D, S, R}(
                                        equation, tableau, params, solver, iguess, cache)
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,TT,D,S,R},
                                        params::ParametersVSPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[3*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[3*(D*(i-1)+k-1)+2]
            cache.Vi[i][k] = x[3*(D*(i-1)+k-1)+3]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]
        params.f_f(tpᵢ, cache.Qi[i], cache.Vi[i], cache.Fi[i])
        params.f_p(tpᵢ, cache.Qi[i], cache.Vi[i], cache.Φi[i])
        cache.Φi[i] .-= cache.Pi[i]
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.tab.λ.c[i]
        params.f_u(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.f_g(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.f_ϕ(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end

    if length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[3*D*S+3*D*R+k]
        end
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.q̃ .+= params.Δt .* params.tab.q.b[i] .* cache.Vi[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[i] .* cache.Fi[i]
    end
    for i in 1:R
        cache.q̃ .+= params.Δt .* params.tab.q.β[i] .* cache.Up[i]
        cache.p̃ .+= params.Δt .* params.tab.p.β[i] .* cache.Gp[i]
    end

    # compute ϕ(q,p)
    tλᵢ = params.t + params.Δt
    params.f_ϕ(tλᵢ, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
@generated function function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersVSPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    cache = IntegratorCacheSPARK{ST,TT,D,S,R}()

    quote
        compute_stages!(y, $cache, params)

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
        for i in 1:S
            for k in 1:D
                b[3*(D*(i-1)+k-1)+1] = - $cache.Yi[i][k]
                b[3*(D*(i-1)+k-1)+2] = - $cache.Zi[i][k]
                b[3*(D*(i-1)+k-1)+3] = - $cache.Φi[i][k]
                for j in 1:S
                    b[3*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * $cache.Vi[j][k]
                    b[3*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[3*(D*(i-1)+k-1)+1] += params.tab.q.α[i,j] * $cache.Up[j][k]
                    b[3*(D*(i-1)+k-1)+2] += params.tab.p.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [(Y-AV-AU), (Z-AF-AG), ωΦ]
        for i in 1:R
            for k in 1:D
                b[3*D*S+3*(D*(i-1)+k-1)+1] = - $cache.Yp[i][k]
                b[3*D*S+3*(D*(i-1)+k-1)+2] = - $cache.Zp[i][k]
                for j in 1:S
                    b[3*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.a[i,j] * $cache.Vi[j][k]
                    b[3*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[3*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.α[i,j] * $cache.Up[j][k]
                    b[3*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.α[i,j] * $cache.Gp[j][k]
                end
            end
        end
        for i in 1:R-P
            for k in 1:D
                b[3*D*S+3*(D*(i-1)+k-1)+3] = 0
                for j in 1:R
                    b[3*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.ω[i,j] * $cache.Φp[j][k]
                end
                b[3*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.ω[i,R+1] * $cache.ϕ̃[k]
            end
        end

        # compute b = d_λ ⋅ Λ
        for i in R-P+1:R
            for k in 1:D
                b[3*D*S+3*(D*(R-1)+k-1)+3] = 0
                for j in 1:R
                    b[3*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.δ[j] * $cache.Λp[j][k]
                end
            end
        end

        if length(params.tab.d) > 0
            for i in 1:S
                for k in 1:D
                    b[3*(D*(i-1)+k-1)+3] -= $cache.μ[k] * params.tab.d[i]
                end
            end

            for k in 1:D
                b[3*D*S+3*D*R+k] = 0
                for i in 1:S
                    b[3*D*S+3*D*R+k] -= $cache.Vi[i][k] * params.tab.d[i]
                end
            end
        end
    end
end


text:
I = ( pi/4 + atan(3.0) ) / 2
for endpt in [neither, both, left, right]
    x, w = legendre(40, endpt)
    Q = sum(w .* 1./(1+(2*x-1).^2))
    @test_approx_eq I Q
end


text:
#=
Problem 3
The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143 ?
=#

import Primes.factor

target = 600851475143
factors = collect(keys(factor(target)))

maxFact = maximum(factors)

println("The greatest prime factor is: $(maxFact)")

max_stars_repo_name: mebrunet/project-euler

max_issues_repo_name: mebrunet/project-euler

max_forks_repo_name: mebrunet/project-euler


text:
# Test pair-wise correlation and covariance functions
x, y = rand(1000), rand(1000)

@test nancov(x,y) ≈ cov(x,y)
@test nancov(x,y, corrected=false) ≈ cov(x,y, corrected=false)
@test nancor(x,y) ≈ cor(x,y)

# Test correlation and covariance functions as applied to matrices
X = rand(100,10)

@test nancov(X) ≈ cov(X)
@test nancov(X, dims=1) ≈ cov(X, dims=1)
@test nancov(X, dims=2) ≈ cov(X, dims=2)
@test nancov(X, corrected=false) ≈ cov(X, corrected=false)

@test nancor(X) ≈ cor(X)
@test nancor(X, dims=1) ≈ cor(X, dims=1)
@test nancor(X, dims=2) ≈ cor(X, dims=2)

# Test NaN-ignoring
xn = [1,2,3,NaN]
xnn = [1,2,3]
@test nancov(xn,xn) ≈ cov(xnn,xnn) ≈ 1
@test nancor(xn,xn) ≈ cor(xnn,xnn) ≈ 1


text:
"""

Returns parameters of a two phase hyper-exponential fitting a mean and an SCV.

"""
function hyper_exp_init(mean_desired::Float64, scv_desired::Float64)::PHDist
    scv_desired < 1.0 && error("SCV must be greater than 1")
    μ1 = 1/(scv_desired+1) #mean parameter
    p = (scv_desired-1)/(scv_desired+1+2/(μ1^2)-4/μ1)
    μ2 = (1-p)/(1-p/μ1) #mean parameter
    α = zeros(2)'
    α[1] = p
    α[2] = 1-p

    T = zeros(2,2)
    T[1,1] = -1/μ1
    T[2,2] = -1/μ2
    return PHDist(α, (1/mean_desired)*T^(-1))
end

"""

Returns parameters of a hypo-exponential (generalized erlang) dist which is a sum of n exponentials with the last one different

"""
function hypo_exp_init(mean::Float64,scv::Float64)::PHDist

    scv ≥ 1.0 && error("SCV must be less than 1")

    n = Int(ceil(1/scv))

    ν1 = n/(1+sqrt((n-1)*(n*scv-1)))
    # ν2 = ν1*(n-1)/(ν1-1)
    ν2 = -(n-1)/(1-ν1)

    α = zeros(n)'
    α[1] = 1
    T = zeros(n,n)
    T[1,1] = -ν1
    T[1,2] = ν1

    for i = 2:(n-1)
        T[i,i] = -ν2
        T[i,i+1] = ν2
    end

    T[n,n] = -ν2

    return PHDist(α, (1/mean)*T)
end


text:
mutable struct outOutput
    outid
    keepid
    outval
    nout
    method
end

mutable struct binomciOutput
    p_hat::Number
    confint::Array
    n::Int
end

mutable struct testOutput
    method
    h0
    n
    df
    estimate
    se
    statistic
    crit
    ci
    p
end

testOutput()=testOutput(nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing)


text:
using TightlyBound
using Test

#include("crystal_testing.jl")
#include("test_ewald.jl")
#include("sym_test.jl")
#include("test_laguerre.jl")
#include("test_U.jl")
#include("test_scf.jl")
#include("test_forces.jl")
#include("test_examples.jl")
#include("test_proto.jl")
#include("test_makescf.jl")
#include("test_atomicproj.jl")

Nothing


text:
@testset "assemble_flux_matrix.jl" begin

end


text:
"""
    update()

Lorem ipsum dolor sit amet.
"""
function update(args::Dict)
  update()
end

function update()
  Pkg.update()
end


text:
using Test
using Klimakoffer
using Printf

EXAMPLES_DIR = joinpath(pathof(Klimakoffer) |> dirname |> dirname, "examples")

@time @testset "Klimakoffer" begin
  test_file = "equilibrium_temperature_1950.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(GlobTemp, 14.484963368770806, atol=1e-12)
  end

  test_file = "equilibrium_temperature_1950_150x75.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(GlobTemp, 14.967892384136475, atol=1e-10)
  end

  test_file = "land_coverage_from_scaled_map.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(landcover, 23.341346153846153, atol=1.5) # aktuell ist es der Wert der originalen 128x65 Karte aus dem Klimakoffer
  end

  test_file = "outline_from_world.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(res, 0, atol=0)
  end

  test_file = "image_to_world.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(result, 0, atol=0.5)
  end

  test_file = "clear_map.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(res, 0, atol=0)
  end

  test_file = "upscale_world.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(res, 0, atol=0)
  end

  test_file = "upscale_albedo.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(res, 0, atol=0)
  end

  test_file = "equilibrium_temperature_albedo.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(GlobTemp, 15.007699256045628, atol=1e-12)
  end

  test_file = "equilibrium_temperature_co2.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(GlobTemp, 14.490442457180537, atol=1e-12)
  end

  test_file = "transient_temperature_co2.txt"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(sol.mean_temperature_yearly[end], 15.124553126197249, atol=1e-12)
  end

  test_file = "equilibrium_temperature_sea_ice_extent.jl"
  @testset "$test_file" begin
    println("")
    println("Running ",test_file)
    println("")
    @test_nowarn include(joinpath(EXAMPLES_DIR, test_file))

    @test isapprox(GlobTemp, 13.552715282159404, atol=1e-12)
  end

  @testset "Printing types to the REPL" begin
    mesh = Mesh()
    @test_nowarn show(stdout, mesh)
    println(stdout)

    model = Model(mesh, 48)
    @test_nowarn show(stdout, model)
    println(stdout)

    discretization = Discretization(mesh, model, 48)
    @test_nowarn show(stdout, discretization)
    println(stdout)
  end
end


text:
module NewtonsMethod

using ForwardDiff
using LinearAlgebra

greet() = print("Hello World!")

include("newtonroot.jl")

export newtonroot
end


text:
@inline symmetric(X::AbstractMatrix) = Symmetric(X)
@inline symmetric(X::Diagonal) = X

diag_Xt_invA_X(A::Cholesky, X::AbstractVecOrMat) = AbstractGPs.diag_At_A(A.U' \ X)

Xt_invA_X(A::Cholesky, x::AbstractVector) = sum(abs2, A.U' \ x)

function Xt_invA_X(A::Cholesky, X::AbstractMatrix)
    V = A.U' \ X
    return Symmetric(V'V)
end

function diag_At_B(A::AbstractVecOrMat, B::AbstractVecOrMat)
    @assert size(A) == size(B)
    return vec(sum(A .* B; dims=1))
end


text:
using OmegaCore
using Distributions
using Test
using OmegaTest

@testset "models" begin
  μ = 1 ~ Uniform(1.0, 2.0)
  x = 2 ~ Normalₚ(μ, 1)
  @test isinferred(randsample, x)

  μ = 1 ~ Normalₚ(0, 1)
  y = 2 ~ Normalₚ(μ, 1)
  @test isinferred(randsample, y)
end


text:
using CLIMAParameters.Planet: R_d, cv_d, T_0, e_int_v0, e_int_i0

"""
    linearized_air_pressure(ρ, ρe_tot, ρe_pot, ρq_tot=0, ρq_liq=0, ρq_ice=0)

The air pressure, linearized around a dry rest state, from the equation of state
(ideal gas law) where:

 - `ρ` (moist-)air density
 - `ρe_tot` total energy density
 - `ρe_pot` potential energy density
and, optionally,
 - `ρq_tot` total water density
 - `ρq_liq` liquid water density
 - `ρq_ice` ice density
"""
function linearized_air_pressure(
    param_set::AbstractParameterSet,
    ρ::FT,
    ρe_tot::FT,
    ρe_pot::FT,
    ρq_tot::FT = FT(0),
    ρq_liq::FT = FT(0),
    ρq_ice::FT = FT(0),
) where {FT <: Real, PS}
    _R_d::FT = R_d(param_set)
    _cv_d::FT = cv_d(param_set)
    _T_0::FT = T_0(param_set)
    _e_int_v0::FT = e_int_v0(param_set)
    _e_int_i0::FT = e_int_i0(param_set)
    return ρ * _R_d * _T_0 +
           _R_d / _cv_d * (
        ρe_tot - ρe_pot - (ρq_tot - ρq_liq) * _e_int_v0 +
        ρq_ice * (_e_int_i0 + _e_int_v0)
    )
end

@inline function linearized_pressure(
    ::DryModel,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    return linearized_air_pressure(param_set, state.ρ, state.ρe, ρe_pot)
end
@inline function linearized_pressure(
    ::EquilMoist,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(
        param_set,
        state.ρ,
        state.ρe,
        ρe_pot,
        state.moisture.ρq_tot,
    )
end

abstract type AtmosLinearModel <: BalanceLaw end

function vars_state_conservative(lm::AtmosLinearModel, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        turbulence::vars_state_conservative(lm.atmos.turbulence, FT)
        hyperdiffusion::vars_state_conservative(lm.atmos.hyperdiffusion, FT)
        moisture::vars_state_conservative(lm.atmos.moisture, FT)
    end
end
vars_state_gradient(lm::AtmosLinearModel, FT) = @vars()
vars_state_gradient_flux(lm::AtmosLinearModel, FT) = @vars()
vars_state_auxiliary(lm::AtmosLinearModel, FT) =
    vars_state_auxiliary(lm.atmos, FT)
vars_integrals(lm::AtmosLinearModel, FT) = @vars()
vars_reverse_integrals(lm::AtmosLinearModel, FT) = @vars()

function update_auxiliary_state!(
    dg::DGModel,
    lm::AtmosLinearModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    return false
end
function flux_second_order!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
integral_load_auxiliary_state!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
integral_set_auxiliary_state!(lm::AtmosLinearModel, aux::Vars, integ::Vars) =
    nothing
reverse_integral_load_auxiliary_state!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
reverse_integral_set_auxiliary_state!(
    lm::AtmosLinearModel,
    aux::Vars,
    integ::Vars,
) = nothing
flux_second_order!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
) = nothing
function wavespeed(lm::AtmosLinearModel, nM, state::Vars, aux::Vars, t::Real)
    ref = aux.ref_state
    return soundspeed_air(lm.atmos.param_set, ref.T)
end

function boundary_state!(
    nf::NumericalFluxFirstOrder,
    atmoslm::AtmosLinearModel,
    args...,
)
    atmos_boundary_state!(nf, AtmosBC(), atmoslm, args...)
end
function boundary_state!(
    nf::NumericalFluxSecondOrder,
    atmoslm::AtmosLinearModel,
    args...,
)
    nothing
end
init_state_auxiliary!(lm::AtmosLinearModel, aux::Vars, geom::LocalGeometry) =
    nothing
init_state_conservative!(
    lm::AtmosLinearModel,
    state::Vars,
    aux::Vars,
    coords,
    t,
) = nothing

struct AtmosAcousticLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("AtmosAcousticLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end

function flux_first_order!(
    lm::AtmosAcousticLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)

    flux.ρ = state.ρu
    pL = linearized_pressure(
        lm.atmos.moisture,
        lm.atmos.param_set,
        lm.atmos.orientation,
        state,
        aux,
    )
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ - e_pot) * state.ρu
    nothing
end
source!(::AtmosAcousticLinearModel, _...) = nothing

struct AtmosAcousticGravityLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticGravityLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("AtmosAcousticGravityLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end
function flux_first_order!(
    lm::AtmosAcousticGravityLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)

    flux.ρ = state.ρu
    pL = linearized_pressure(
        lm.atmos.moisture,
        lm.atmos.param_set,
        lm.atmos.orientation,
        state,
        aux,
    )
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ) * state.ρu
    nothing
end
function source!(
    lm::AtmosAcousticGravityLinearModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    if direction isa VerticalDirection || direction isa EveryDirection
        ∇Φ = ∇gravitational_potential(lm.atmos.orientation, aux)
        source.ρu -= state.ρ * ∇Φ
    end
    nothing
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA