text:
function randmech()
    rand_tree_mechanism(Float64,
        QuaternionFloating{Float64},
        [Revolute{Float64} for i = 1 : 5]...,
        [Fixed{Float64} for i = 1 : 5]...,
        [Prismatic{Float64} for i = 1 : 5]...,
        [Planar{Float64} for i = 1 : 5]...,
        [SPQuatFloating{Float64} for i = 1:2]...,
        [SinCosRevolute{Float64} for i = 1:2]...
    )
end

@testset "mechanism algorithms" begin
    @testset "show" begin
        Random.seed!(17)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        mechanism_with_loops = deepcopy(mechanism)
        for i = 1 : 5
            pred = rand(bodies(mechanism_with_loops))
            succ = rand(bodies(mechanism_with_loops))
            joint = Joint("non-tree-$i", Fixed{Float64}())
            attach!(mechanism_with_loops, pred, succ, joint)
        end

        show(devnull, mechanism_with_loops)
        for joint in joints(mechanism_with_loops)
            show(devnull, joint)
            show(IOContext(devnull, :compact => true), joint)
        end
        for body in bodies(mechanism_with_loops)
            show(devnull, body)
            show(IOContext(devnull, :compact => true), body)
        end
        show(devnull, x)
    end

    @testset "supports" begin
        Random.seed!(25)
        mechanism = randmech()
        mc_mechanism = maximal_coordinates(mechanism)
        for m in [mechanism, mc_mechanism]
            state = MechanismState(m)
            for body in bodies(m)
                body_ancestors = RigidBodyDynamics.Graphs.ancestors(body, m.tree)
                for joint in tree_joints(m)
                    @test RigidBodyDynamics.supports(joint, body, state) == (successor(joint, m) ∈ body_ancestors)
                end
                for joint in non_tree_joints(m)
                    @test !RigidBodyDynamics.supports(joint, body, state)
                end
            end
        end
    end

    @testset "basic stuff" begin
        Random.seed!(18)
        mechanism = randmech()
        for joint in joints(mechanism)
            @test eltype(joint_type(joint)) == Float64
        end

        x = MechanismState(mechanism)
        rand!(x)
        q = vcat([configuration(x, joint) for joint in tree_joints(mechanism)]...)
        v = vcat([velocity(x, joint) for joint in tree_joints(mechanism)]...)

        @test q == configuration(x)
        @test v == velocity(x)

        zero_configuration!(x)
        set_configuration!(x, q)
        @test q == configuration(x)

        zero_velocity!(x)
        set_velocity!(x, v)
        @test v == velocity(x)

        qcopy = copy(configuration(x))
        zero_configuration!(x)
        for joint in joints(mechanism)
            set_configuration!(x, joint, qcopy[configuration_range(x, joint)])
        end
        @test q == configuration(x)

        vcopy = copy(velocity(x))
        zero_velocity!(x)
        for joint in joints(mechanism)
            set_velocity!(x, joint, vcopy[velocity_range(x, joint)])
        end
        @test v == velocity(x)

        zero!(x)
        copyto!(x, [q; v])

        @test q == configuration(x)
        @test v == velocity(x)

        q2 = rand(num_positions(mechanism))
        v2 = rand(num_velocities(mechanism))
        q2copy = deepcopy(q2)
        v2copy = deepcopy(v2)
        x2 = MechanismState(mechanism, q2, v2)
        @test parent(configuration(x2)) === q2
        @test parent(velocity(x2)) === v2
        @test all(configuration(x2) .== q2copy)
        @test all(velocity(x2) .== v2copy)

        @test @inferred(num_positions(mechanism)) == num_positions(x)
        @test @inferred(num_velocities(mechanism)) == num_velocities(x)

        for joint in tree_joints(mechanism)
            for i in configuration_range(x, joint)
                @test RigidBodyDynamics.configuration_index_to_joint_id(x, i) == JointID(joint)
            end
            for i in velocity_range(x, joint)
                @test RigidBodyDynamics.velocity_index_to_joint_id(x, i) == JointID(joint)
            end
        end

        rand!(x)
        q = configuration(x)
        @test q != qcopy
        for joint in tree_joints(mechanism)
            q[joint] .= qcopy[joint]
        end
        @test q == qcopy
    end

    @testset "copyto! / Vector" begin
        Random.seed!(19)
        mechanism = randmech()
        x1 = MechanismState(mechanism)
        rand!(x1)
        @test Vector(x1) == [configuration(x1); velocity(x1); additional_state(x1)]

        x2 = MechanismState(mechanism)
        rand!(x2)
        copyto!(x1, x2)
        @test Vector(x1) == Vector(x2)

        x3 = MechanismState(mechanism)
        copyto!(x3, Vector(x2))
        @test Vector(x3) == Vector(x2)

        @test Vector{Float32}(x1) isa Vector{Float32}
        @test Vector{Float32}(x1) ≈ Vector(x1) atol=1e-6
        @test Array(x1) == Array{Float64}(x1) == convert(Array, x1) == convert(Array{Float64}, x1)
        @test Vector(x1) == Vector{Float64}(x1) == convert(Vector, x1) == convert(Vector{Float64}, x1)
    end

    @testset "q̇ <-> v" begin
        Random.seed!(20)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        q = configuration(x)
        q̇ = configuration_derivative(x)
        v = velocity(x)
        for joint in joints(mechanism)
            qjoint = configuration(x, joint)
            q̇joint = q̇[configuration_range(x, joint)]
            vjoint = velocity(x, joint)
            vjoint_from_q̇ = similar(vjoint)
            configuration_derivative_to_velocity!(vjoint_from_q̇, joint, qjoint, q̇joint)
            @test isapprox(vjoint, vjoint_from_q̇; atol = 1e-12)
        end

        Jq̇_to_v = RigidBodyDynamics.configuration_derivative_to_velocity_jacobian(x)
        Jv_to_q̇ = RigidBodyDynamics.velocity_to_configuration_derivative_jacobian(x)
        for i in 1:10
            rand!(x)
            q = configuration(x)
            q̇ = configuration_derivative(x)
            v = velocity(x)

            RigidBodyDynamics.configuration_derivative_to_velocity_jacobian!(Jq̇_to_v, x)
            RigidBodyDynamics.velocity_to_configuration_derivative_jacobian!(Jv_to_q̇, x)
            @test Jq̇_to_v * q̇ ≈ v
            @test Jv_to_q̇ * v ≈ q̇

            @test @allocated(RigidBodyDynamics.configuration_derivative_to_velocity_jacobian!(Jq̇_to_v, x)) == 0
            @test @allocated(RigidBodyDynamics.velocity_to_configuration_derivative_jacobian!(Jv_to_q̇, x)) == 0

            for joint in joints(mechanism)
                qrange = configuration_range(x, joint)
                vrange = velocity_range(x, joint)
                qj = q[qrange]
                q̇j = q̇[qrange]
                vj = v[vrange]

                Jv_to_q̇_j = RigidBodyDynamics.velocity_to_configuration_derivative_jacobian(joint, qj)
                Jq̇_to_v_j = RigidBodyDynamics.configuration_derivative_to_velocity_jacobian(joint, qj)

                if num_velocities(joint) > 0
                    @test Jv_to_q̇_j * vj ≈ q̇j
                    @test Jq̇_to_v_j * q̇j ≈ vj
                else
                    @test size(Jv_to_q̇_j) == (0, 0)
                    @test size(Jq̇_to_v_j) == (0, 0)
                end
            end
        end
    end

    @testset "set_configuration! / set_velocity!" begin
        Random.seed!(21)
        mechanism = randmech()
        x = MechanismState(mechanism)
        for joint in joints(mechanism)
            qjoint = rand(num_positions(joint))
            set_configuration!(x, joint, qjoint)
            @test configuration(x, joint) == qjoint
            @test configuration(x, joint) !== qjoint

            vjoint = rand(num_velocities(joint))
            set_velocity!(x, joint, vjoint)
            @test velocity(x, joint) == vjoint
            @test velocity(x, joint) !== vjoint

            if joint_type(joint) isa QuaternionFloating
                tf = rand(Transform3D{Float64}, frame_after(joint), frame_before(joint))
                set_configuration!(x, joint, tf)
                @test RigidBodyDynamics.joint_transform(joint, configuration(x, joint)) ≈ tf atol = 1e-12

                # TODO: the frame stuff is kind of awkward here.
                twist = RigidBodyDynamics.joint_twist(joint, configuration(x, joint), velocity(x, joint))
                twist = rand(Twist{Float64}, frame_after(joint), frame_before(joint), frame_after(joint))
                set_velocity!(x, joint, twist)
                twist_back = RigidBodyDynamics.joint_twist(joint, configuration(x, joint), velocity(x, joint))
                @test twist_back.angular ≈ twist.angular atol = 1e-12
                @test twist_back.linear ≈ twist.linear atol = 1e-12
            end
            if joint_type(joint) isa Revolute || joint_type(joint) isa Prismatic
                qj = rand()
                set_configuration!(x, joint, qj)
                @test configuration(x, joint)[1] == qj

                vj = rand()
                set_velocity!(x, joint, vj)
                @test velocity(x, joint)[1] == vj
            end
            if joint_type(joint) isa QuaternionSpherical
                quat = rand(UnitQuaternion{Float64})
                set_configuration!(x, joint, quat)
                tf = RigidBodyDynamics.joint_transform(joint, configuration(x, joint))
                @test UnitQuaternion(rotation(tf)) ≈ quat atol = 1e-12
            end
            if joint_type(joint) isa SinCosRevolute
                qj = rand()
                set_configuration!(x, joint, qj)
                @test SVector(sincos(qj)) == configuration(x, joint)

                vj = rand()
                set_velocity!(x, joint, vj)
                @test velocity(x, joint)[1] == vj
            end
        end
    end

    @testset "normalize_configuration!" begin
        Random.seed!(22)
        mechanism = randmech()
        let x = MechanismState(mechanism) # required to achieve zero allocations
            configuration(x) .= 1
            for joint in joints(mechanism)
                qjoint = configuration(x, joint)
                requires_normalization = num_positions(joint) != num_velocities(joint) # TODO: not quite the right thing to check for
                @test requires_normalization != RigidBodyDynamics.is_configuration_normalized(joint, qjoint)
            end
            normalize_configuration!(x)
            for joint in joints(mechanism)
                @test RigidBodyDynamics.is_configuration_normalized(joint, configuration(x, joint))
            end
            allocs = @allocated normalize_configuration!(x)
            @test allocs == 0
        end
    end

    @testset "joint_torque! / motion_subspace" begin
        Random.seed!(23)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        for joint in tree_joints(mechanism)
            body = successor(joint, mechanism)
            qjoint = configuration(x, joint)
            wrench = rand(Wrench{Float64}, frame_after(joint))
            τ = Vector{Float64}(undef, num_velocities(joint))
            RigidBodyDynamics.joint_torque!(τ, joint, qjoint, wrench)
            S = motion_subspace(joint, configuration(x, joint))
            @test isapprox(τ, torque(S, wrench))
        end
    end

    @testset "isfloating" begin
        Random.seed!(24)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        for joint in joints(mechanism)
            num_positions(joint) == 0 && continue # https://github.com/JuliaLang/julia/issues/26578
            S = motion_subspace(joint, configuration(x, joint))
            @test isfloating(joint) == (rank(Array(S)) == 6)
            @test isfloating(joint) == (num_constraints(joint) == 0)
        end
    end

    @testset "geometric_jacobian / relative_twist" begin
        Random.seed!(25)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        frame = CartesianFrame3D()
        for i = 1 : 100
            bs = Set(bodies(mechanism))
            body = rand([bs...])
            delete!(bs, body)
            base = rand([bs...])
            p = path(mechanism, base, body)
            v = velocity(x)

            J = geometric_jacobian(x, p)
            T = relative_twist(x, body, base)
            @test isapprox(Twist(J, v), T; atol = 1e-12)

            J1 = GeometricJacobian(J.body, J.base, J.frame, similar(angular(J)), similar(linear(J)))
            geometric_jacobian!(J1, x, p)
            @test isapprox(Twist(J1, v), T; atol = 1e-12)

            H = rand(Transform3D, root_frame(mechanism), frame)
            J2 = GeometricJacobian(J.body, J.base, frame, similar(angular(J)), similar(linear(J)))
            if num_velocities(p) > 0
                @test_throws ArgumentError geometric_jacobian!(J, x, p, H)
            end
            geometric_jacobian!(J2, x, p, H)
            @test isapprox(Twist(J2, v), transform(T, H); atol = 1e-12)

            J3 = GeometricJacobian(J.body, J.base, default_frame(body), similar(angular(J)), similar(linear(J)))
            geometric_jacobian!(J3, x, p)
            @test isapprox(Twist(J3, v), transform(x, T, default_frame(body)); atol = 1e-12)
        end
    end

    @testset "point jacobian" begin
        Random.seed!(26)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        @testset "point expressed in body frame" begin
            for i = 1 : 100
                bs = Set(bodies(mechanism))
                body = rand([bs...])
                delete!(bs, body)
                base = rand([bs...])
                p = path(mechanism, base, body)
                point = Point3D(default_frame(body), rand(SVector{3, Float64}))
                J_point = point_jacobian(x, p, point)

                # Check agreement with GeometricJacobian -> Twist -> point_velocity
                J = geometric_jacobian(x, p)
                T = Twist(J, velocity(x))
                point_in_world = transform(x, point, root_frame(mechanism))
                point_velocity_expected = point_velocity(T, point_in_world)
                @test point_velocity_expected ≈ transform(x, point_velocity(J_point, velocity(x)), root_frame(mechanism))

                # Test that point_velocity give us what Jp * v does
                @test transform(x, point_velocity(J_point, velocity(x)), point.frame).v ≈ Array(J_point) * velocity(x)

                # Test that in-place updates work too
                rand!(x)
                if point.frame != default_frame(base)
                    @test_throws ArgumentError point_jacobian!(J_point, x, p, transform(x, point, default_frame(base)))
                end
                point_jacobian!(J_point, x, p, point)
                @test(@ballocated(point_jacobian!($J_point, $x, $p, $point)) == 0)
                J = geometric_jacobian(x, p)
                T = Twist(J, velocity(x))
                point_in_world = transform(x, point, root_frame(mechanism))
                @test point_velocity(T, point_in_world) ≈ transform(x, point_velocity(J_point, velocity(x)), root_frame(mechanism))

                # Test Jᵀ * f
                f = FreeVector3D(CartesianFrame3D(), rand(), rand(), rand())
                τ = similar(velocity(x))
                @test_throws ArgumentError mul!(τ, transpose(J_point), f)
                f = FreeVector3D(J_point.frame, rand(), rand(), rand())
                mul!(τ, transpose(J_point), f)
                @test τ == transpose(J_point.J) * f.v
                @test τ == transpose(J_point) * f
            end
        end

        @testset "point expressed in world frame" begin
            Random.seed!(27)
            for i = 1 : 10
                bs = Set(bodies(mechanism))
                body = rand([bs...])
                delete!(bs, body)
                base = rand([bs...])
                p = path(mechanism, base, body)
                point = Point3D(root_frame(mechanism), rand(SVector{3, Float64}))
                J_point = point_jacobian(x, p, point)

                # Check agreement with GeometricJacobian -> Twist -> point_velocity
                J = geometric_jacobian(x, p)
                T = Twist(J, velocity(x))
                @test point_velocity(T, point) ≈ point_velocity(J_point, velocity(x))

                # Test that point_velocity give us what Jp * v does
                @test point_velocity(J_point, velocity(x)).v ≈ Array(J_point) * velocity(x)

                point_jacobian!(J_point, x, p, point)
            end
        end

    end

    @testset "motion_subspace / constraint_wrench_subspace" begin
        Random.seed!(28)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        for joint in tree_joints(mechanism)
            body = successor(joint, mechanism)
            qjoint = configuration(x, joint)
            S = motion_subspace(joint, configuration(x, joint))
            tf = joint_transform(joint, qjoint)
            T = constraint_wrench_subspace(joint, tf)
            if 0 < num_constraints(joint) < 6
                @test isapprox(angular(T)' * angular(S) + linear(T)' * linear(S), zeros(num_constraints(joint), num_velocities(joint)); atol = 1e-12)
            elseif num_constraints(joint) == 0
                @test size(T) == (6, 0)
            else
                @test size(S) == (6, 0)
            end
        end
    end

    # TODO: good test, but currently don't want to support joints with variable wrench subspaces:
    # @testset "constraint_bias" begin
    #     for joint in joints(mechanism)
    #         qjoint = configuration(x, joint)
    #         vjoint = velocity(x, joint)
    #         q̇joint = similar(qjoint)
    #         velocity_to_configuration_derivative!(joint, q̇joint, qjoint, vjoint)
    #         qjoint_autodiff = ForwardDiff.Dual.(qjoint, q̇joint)
    #         TAutodiff = constraint_wrench_subspace(joint, qjoint_autodiff)#::RigidBodyDynamics.WrenchSubspace{eltype(qjoint_autodiff)} # TODO
    #         ang = map(x -> ForwardDiff.partials(x, 1), angular(TAutodiff))
    #         lin = map(x -> ForwardDiff.partials(x, 1), linear(TAutodiff))
    #         Ṫ = WrenchMatrix(frame_after(joint), ang, lin)
    #         joint_twist = transform(x, relative_twist(x, frame_after(joint), frame_before(joint)), frame_after(joint))
    #         bias = fill(NaN, 6 - num_velocities(joint))
    #         constraint_bias!(joint, bias, joint_twist)
    #         @test isapprox(angular(Ṫ)' * angular(joint_twist) + linear(Ṫ)' * linear(joint_twist), bias; atol = 1e-14)
    #     end
    # end

    @testset "relative_acceleration" begin
        Random.seed!(29)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        result = DynamicsResult(mechanism)
        for body in bodies(mechanism)
            for base in bodies(mechanism)
                v̇ = similar(velocity(x))
                rand!(v̇)
                spatial_accelerations!(result.accelerations, x, v̇)
                Ṫ = relative_acceleration(result, body, base)
                q = configuration(x)
                v = velocity(x)
                q̇ = configuration_derivative(x)
                q_autodiff = ForwardDiff.Dual.(q, q̇)
                v_autodiff = ForwardDiff.Dual.(v, v̇)
                x_autodiff = MechanismState{eltype(q_autodiff)}(mechanism)
                set_configuration!(x_autodiff, q_autodiff)
                set_velocity!(x_autodiff, v_autodiff)
                twist_autodiff = relative_twist(x_autodiff, body, base)
                accel_vec = [ForwardDiff.partials(x, 1)::Float64 for x in (SVector(twist_autodiff))]
                @test isapprox(SVector(Ṫ), accel_vec; atol = 1e-12)

                root = root_body(mechanism)
                f = default_frame(body)
                Ṫbody = transform(x, relative_acceleration(result, body, root), f)
                Ṫbase = transform(x, relative_acceleration(result, base, root), f)
                @test isapprox(transform(x, -Ṫbase + Ṫbody, Ṫ.frame), Ṫ; atol = 1e-12)
            end
        end
    end

    @testset "motion subspace / twist wrt world" begin
        Random.seed!(30)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        for joint in tree_joints(mechanism)
            body = successor(joint, mechanism)
            parent_body = predecessor(joint, mechanism)
            toroot = transform_to_root(x, body)
            S = transform(motion_subspace(joint, configuration(x, joint)), toroot)
            @test isapprox(relative_twist(x, body, parent_body), Twist(S, velocity(x, joint)); atol = 1e-12)
        end
    end

    @testset "composite rigid body inertias" begin
        Random.seed!(31)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        for joint in tree_joints(mechanism)
            body = successor(joint, mechanism)
            crb = crb_inertia(x, body)
            stack = [body]
            subtree = typeof(body)[]
            while !isempty(stack)
                b = pop!(stack)
                push!(subtree, b)
                for joint in out_joints(b, mechanism)
                    push!(stack, successor(joint, mechanism))
                end
            end
            @test isapprox(sum((b::RigidBody) -> spatial_inertia(x, b), subtree), crb; atol = 1e-12)
        end
    end

    @testset "momentum_matrix / summing momenta" begin
        Random.seed!(32)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        A = momentum_matrix(x)
        Amat = Array(A)
        for joint in tree_joints(mechanism)
            if num_velocities(joint) > 0 # TODO: Base.mapslices can't handle matrices with zero columns
                body = successor(joint, mechanism)
                Ajoint = Amat[:, velocity_range(x, joint)]
                S = transform(motion_subspace(joint, configuration(x, joint)), transform_to_root(x, body))
                @test isapprox(Array(crb_inertia(x, body) * S), Ajoint; atol = 1e-12)
            end
        end

        v = velocity(x)
        hsum = sum(b -> spatial_inertia(x, b) * twist_wrt_world(x, b), non_root_bodies(mechanism))
        @test isapprox(Momentum(A, v), hsum; atol = 1e-12)

        A1 = MomentumMatrix(A.frame, similar(angular(A)), similar(linear(A)))
        momentum_matrix!(A1, x)
        @test isapprox(Momentum(A1, v), hsum; atol = 1e-12)

        frame = CartesianFrame3D()
        A2 = MomentumMatrix(frame, similar(angular(A)), similar(linear(A)))
        H = rand(Transform3D, root_frame(mechanism), frame)
        @test_throws ArgumentError momentum_matrix!(A, x, H)
        momentum_matrix!(A2, x, H)
        @test isapprox(Momentum(A2, v), transform(hsum, H); atol = 1e-12)

        body = rand(collect(bodies(mechanism)))
        A3 = MomentumMatrix(default_frame(body), similar(angular(A)), similar(linear(A)))
        momentum_matrix!(A3, x)
        @test isapprox(Momentum(A3, v), transform(x, hsum, default_frame(body)); atol = 1e-12)
    end

    @testset "mass matrix / kinetic energy" begin
        Random.seed!(33)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        Ek = kinetic_energy(x)
        M = mass_matrix(x)
        v = velocity(x)
        @test isapprox(1/2 * dot(v, M * v), Ek; atol = 1e-12)

        q = configuration(x)
        cache = StateCache(mechanism)
        kinetic_energy_fun = function (v)
            local x = cache[eltype(v)]
            set_configuration!(x, q)
            set_velocity!(x, v)
            kinetic_energy(x)
        end

        # FIXME: excessive compilation time
        # M2 = similar(M.data)
        # NOTE: chunk size 1 necessary after updating to ForwardDiff 0.2 because creating a MechanismState with a max size Dual takes forever...
        # see https://github.com/JuliaDiff/ForwardDiff.jl/issues/266
        # M2 = ForwardDiff.hessian!(M2, kinetic_energy_fun, v, ForwardDiff.HessianConfig(kinetic_energy_fun, v, ForwardDiff.Chunk{1}()))
        # @test isapprox(M2, M; atol = 1e-12)
    end

    @testset "spatial_inertia!" begin
        Random.seed!(34)
        mechanism = randmech()
        body = rand(collect(non_root_bodies(mechanism)))
        newinertia = rand(SpatialInertia{eltype(mechanism)}, spatial_inertia(body).frame)
        spatial_inertia!(body, newinertia)
        @assert spatial_inertia(body) == newinertia
    end

    @testset "inverse dynamics / acceleration term" begin
        Random.seed!(35)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        M = mass_matrix(x)
        function v̇_to_τ(v̇)
            inverse_dynamics(x, v̇)
        end
        M2 = similar(M.data)
        v̇ = similar(velocity(x))
        v̇ .= 0
        ForwardDiff.jacobian!(M2, v̇_to_τ, v̇)
        @test isapprox(M2, M; atol = 1e-12)
    end

    @testset "inverse dynamics / Coriolis term" begin
        Random.seed!(36)
        mechanism = rand_tree_mechanism(Float64, [[Revolute{Float64} for i = 1 : 10]; [Prismatic{Float64} for i = 1 : 10]]...) # skew symmetry property tested later on doesn't hold when q̇ ≠ v
        x = MechanismState{Float64}(mechanism)
        rand!(x)

        cache = StateCache(mechanism)
        function q_to_M(q)
            local x = cache[eltype(q)]
            set_configuration!(x, q)
            zero_velocity!(x)
            vec(mass_matrix(x))
        end
        nv = num_velocities(mechanism)
        nq = num_positions(mechanism)
        dMdq = zeros(nv * nv, nq)
        ForwardDiff.jacobian!(dMdq, q_to_M, configuration(x))
        q̇ = velocity(x)
        Ṁ = reshape(dMdq * q̇, num_velocities(mechanism), num_velocities(mechanism))

        q = configuration(x)
        v̇ = similar(velocity(x))
        v̇ .= 0
        cache = StateCache(mechanism)
        function v_to_c(v)
            local x = cache[eltype(v)]
            set_configuration!(x, q)
            set_velocity!(x, v)
            inverse_dynamics(x, v̇)
        end
        C = similar(Ṁ)
        ForwardDiff.jacobian!(C, v_to_c, q̇)
        C *= 1/2

        skew = Ṁ - 2 * C;
        @test isapprox(skew + skew', zeros(size(skew)); atol = 1e-12)
    end

    @testset "inverse dynamics / gravity term" begin
        Random.seed!(37)
        mechanism = rand_tree_mechanism(Float64, [[Revolute{Float64} for i = 1 : 10]; [Prismatic{Float64} for i = 1 : 10]]...)
        x = MechanismState(mechanism)
        rand!(x)
        v̇ = similar(velocity(x))
        v̇ .= 0
        zero_velocity!(x)
        g = inverse_dynamics(x, v̇)

        cache = StateCache(mechanism)
        function q_to_potential(q)
            local x = cache[eltype(q)]
            set_configuration!(x, q)
            zero_velocity!(x)
            return [gravitational_potential_energy(x)]
        end

        g2 = similar(g')
        ForwardDiff.jacobian!(g2, q_to_potential, configuration(x))
        @test isapprox(g2, g'; atol = 1e-12)
    end

    @testset "momentum matrix" begin
        Random.seed!(38)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        q = configuration(x)
        q̇ = configuration_derivative(x)
        v = velocity(x)
        v̇ = similar(velocity(x))
        rand!(v̇)

        # momentum computed two ways
        @test isapprox(Momentum(momentum_matrix(x), v), momentum(x))

        # rate of change of momentum computed using autodiff:
        q_autodiff = ForwardDiff.Dual.(q, q̇)
        v_autodiff = ForwardDiff.Dual.(v, v̇)
        x_autodiff = MechanismState{eltype(q_autodiff)}(mechanism)
        set_configuration!(x_autodiff, q_autodiff)
        set_velocity!(x_autodiff, v_autodiff)
        A_autodiff = Array(momentum_matrix(x_autodiff))
        A = [ForwardDiff.value(A_autodiff[i, j])::Float64 for i = 1 : size(A_autodiff, 1), j = 1 : size(A_autodiff, 2)]
        Ȧ = [ForwardDiff.partials(A_autodiff[i, j], 1)::Float64 for i = 1 : size(A_autodiff, 1), j = 1 : size(A_autodiff, 2)]
        ḣArray = A * v̇ + Ȧ * v

        # rate of change of momentum computed without autodiff:
        ḣ = Wrench(momentum_matrix(x), v̇) + momentum_rate_bias(x)
        @test isapprox(ḣArray, SVector(ḣ); atol = 1e-12)
    end

    @testset "inverse dynamics / external wrenches" begin
        # test requires a floating mechanism
        Random.seed!(39)
        mechanism = rand_floating_tree_mechanism(Float64, fill(Revolute{Float64}, 10)..., fill(Planar{Float64}, 10)..., fill(SinCosRevolute{Float64}, 5)...)
        x = MechanismState(mechanism)
        rand!(x)
        v̇ = similar(velocity(x))
        rand!(v̇)
        externalwrenches = Dict(BodyID(body) => rand(Wrench{Float64}, root_frame(mechanism)) for body in bodies(mechanism))
        τ = inverse_dynamics(x, v̇, externalwrenches)
        floatingjoint = first(out_joints(root_body(mechanism), mechanism))
        τfloating = τ[velocity_range(x, floatingjoint)]
        floatingjointwrench = Wrench(frame_after(floatingjoint), SVector{3}(τfloating[1 : 3]), SVector{3}(τfloating[4 : 6]))
        floatingjointwrench = transform(x, floatingjointwrench, root_frame(mechanism))
        ḣ = Wrench(momentum_matrix(x), v̇) + momentum_rate_bias(x) # momentum rate of change
        gravitational_force = mass(mechanism) * mechanism.gravitational_acceleration
        com = center_of_mass(x)
        gravitational_wrench = Wrench(gravitational_force.frame, (com × gravitational_force).v, gravitational_force.v)
        total_wrench = floatingjointwrench + gravitational_wrench + sum((b) -> transform(x, externalwrenches[BodyID(b)], root_frame(mechanism)), non_root_bodies(mechanism))
        @test isapprox(total_wrench, ḣ; atol = 1e-10)
    end

    @testset "dynamics / inverse dynamics" begin
        Random.seed!(40)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        external_torques = rand(num_velocities(mechanism))
        externalwrenches = Dict(BodyID(body) => rand(Wrench{Float64}, root_frame(mechanism)) for body in bodies(mechanism))
        result = DynamicsResult(mechanism)
        dynamics!(result, x, external_torques, externalwrenches)
        τ = inverse_dynamics(x, result.v̇, externalwrenches) - external_torques
        @test isapprox(τ, zeros(num_velocities(mechanism)); atol = 1e-10)
    end

    @testset "dynamics_bias / inverse_dynamics" begin
        Random.seed!(41)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        externalwrenches = Dict(BodyID(body) => rand(Wrench{Float64}, root_frame(mechanism)) for body in bodies(mechanism))
        v̇ = similar(velocity(x))
        v̇ .= 0
        τ1 = inverse_dynamics(x, v̇, externalwrenches)
        τ2 = dynamics_bias(x, externalwrenches)
        @test τ1 ≈ τ2
    end

    @testset "dynamics ode method" begin
        Random.seed!(42)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        torques = rand(num_velocities(mechanism))
        externalwrenches = Dict(BodyID(body) => rand(Wrench{Float64}, root_frame(mechanism)) for body in bodies(mechanism))

        result1 = DynamicsResult(mechanism)
        ẋ = similar(Vector(x))
        dynamics!(ẋ, result1, x, Vector(x), torques, externalwrenches)

        result2 = DynamicsResult(mechanism)
        dynamics!(result2, x, torques, externalwrenches)

        @test isapprox([configuration_derivative(x); result2.v̇], ẋ)
    end

    @testset "power flow" begin
        Random.seed!(43)
        mechanism = randmech()
        x = MechanismState(mechanism)
        rand!(x)
        externalwrenches = Dict(BodyID(body) => rand(Wrench{Float64}, root_frame(mechanism)) for body in bodies(mechanism))
        τ = similar(velocity(x))
        rand!(τ)
        result = DynamicsResult(mechanism)
        dynamics!(result, x, τ, externalwrenches)

        q = configuration(x)
        q̇ = configuration_derivative(x)
        v = velocity(x)
        v̇ = result.v̇
        power = τ ⋅ v + sum(body -> externalwrenches[BodyID(body)] ⋅ twist_wrt_world(x, body), non_root_bodies(mechanism))

        q_autodiff = ForwardDiff.Dual.(q, q̇)
        v_autodiff = ForwardDiff.Dual.(v, v̇)
        x_autodiff = MechanismState{eltype(q_autodiff)}(mechanism)
        set_configuration!(x_autodiff, q_autodiff)
        set_velocity!(x_autodiff, v_autodiff)
        energy_autodiff = gravitational_potential_energy(x_autodiff) + kinetic_energy(x_autodiff)
        energy_derivative = ForwardDiff.partials(energy_autodiff)[1]
        @test isapprox(power, energy_derivative, atol = 1e-10)
    end

    @testset "local / global coordinates" begin
        Random.seed!(44)
        mechanism = randmech()
        state = MechanismState(mechanism)
        rand!(state)
        for joint in joints(mechanism)
            # back and forth between local and global
            ϕ = Vector{Float64}(undef, num_velocities(joint))
            ϕ̇ = Vector{Float64}(undef, num_velocities(joint))
            q0 = Vector{Float64}(undef, num_positions(joint))
            q = configuration(state, joint)
            v = velocity(state, joint)
            rand_configuration!(q0, joint)
            local_coordinates!(ϕ, ϕ̇, joint, q0, q, v)
            q_back = Vector{Float64}(undef, num_positions(joint))
            global_coordinates!(q_back, joint, q0, ϕ)
            principal_value!(q_back, joint)

            let expected = copy(q)
                principal_value!(expected, joint)
                @test isapprox(q_back, expected)
            end

            # compare ϕ̇ to autodiff
            q̇ = Vector{Float64}(undef, num_positions(joint))
            velocity_to_configuration_derivative!(q̇, joint, q, v)
            v̇ = rand(num_velocities(joint))
            q_autodiff = ForwardDiff.Dual.(q, q̇)
            v_autodiff = ForwardDiff.Dual.(v, v̇)
            q0_autodiff = ForwardDiff.Dual.(q0, zeros(length(q0)))
            T = eltype(q_autodiff)
            ϕ_autodiff = Vector{T}(undef, length(ϕ))
            ϕ̇_autodiff = Vector{T}(undef, length(ϕ̇))
            local_coordinates!(ϕ_autodiff, ϕ̇_autodiff, joint, q0_autodiff, q_autodiff, v_autodiff)
            ϕ̇_from_autodiff = [ForwardDiff.partials(x)[1] for x in ϕ_autodiff]
            @test isapprox(ϕ̇, ϕ̇_from_autodiff)

            # local coordinates should be zero when q = q0
            # Definition 2.9 in Duindam, "Port-Based Modeling and Control for Efficient Bipedal Walking Robots"
            copyto!(q, q0)
            local_coordinates!(ϕ, ϕ̇, joint, q0, q, v)
            @test isapprox(ϕ, zeros(num_velocities(joint)); atol = 1e-15)
        end
    end

    @testset "configuration_derivative_to_velocity_adjoint!" begin
        Random.seed!(45)
        mechanism = randmech()
        x = MechanismState(mechanism)
        configuration(x) .= rand(num_positions(x)) # needs to work for configuration vectors that do not satisfy the state constraints as well
        fv = similar(velocity(x))
        rand!(fv)
        fq = similar(configuration(x))
        rand!(fq)
        configuration_derivative_to_velocity_adjoint!(fq, x, fv)
        @test dot(fv, velocity(x)) ≈ dot(fq, configuration_derivative(x))
    end

    @testset "joint bounds" begin
        @test @inferred(RigidBodyDynamics.Bounds{Float64}()).lower == -Inf
        @test @inferred(RigidBodyDynamics.Bounds{Float64}()).upper == Inf
        @test @inferred(RigidBodyDynamics.Bounds(-1, 2.0)).lower == -1
        @test @inferred(RigidBodyDynamics.Bounds(-1, 2.0)).upper == 2
        @test isa(RigidBodyDynamics.Bounds(-1, 1.0).lower, Float64)
        @test isa(RigidBodyDynamics.Bounds(-1, 1.0).upper, Float64)
        @test @inferred(clamp(1.5, RigidBodyDynamics.Bounds(-1, 1))) == RigidBodyDynamics.upper(RigidBodyDynamics.Bounds(-1, 1))
        @test @inferred(clamp(-3, RigidBodyDynamics.Bounds(-2, 1))) == RigidBodyDynamics.lower(RigidBodyDynamics.Bounds(-2, 1))
        @test @inferred(clamp(0, RigidBodyDynamics.Bounds(-1, 1))) == 0
        @test @inferred(intersect(RigidBodyDynamics.Bounds(-1, 1), RigidBodyDynamics.Bounds(-0.5, 2.0))) == RigidBodyDynamics.Bounds(-0.5, 1.0)
        @test @inferred(convert(RigidBodyDynamics.Bounds{Float64}, RigidBodyDynamics.Bounds(1, 2))) == RigidBodyDynamics.Bounds{Float64}(1.0, 2.0)
    end

    @testset "issue #330" begin
        Random.seed!(46)
        mechanism = rand_tree_mechanism(Revolute{Float64})
        f330(x::AbstractArray{T}) where {T} = begin
            result = DynamicsResult{T}(mechanism)
            1.
        end
        @test ForwardDiff.hessian(f330, [1.]) == zeros(1, 1)
    end

    @testset "principal_value!" begin
        Random.seed!(47)
        state_orig = MechanismState(randmech())
        Random.rand!(state_orig)
        for joint_k = tree_joints(state_orig.mechanism)
            joint_type_k = joint_type(joint_k)
            if isa(joint_type_k, SPQuatFloating)
                RigidBodyDynamics.set_rotation!(state_orig.q[joint_k], joint_type_k, ModifiedRodriguesParam(UnitQuaternion(-0.5, randn(), randn(), randn())))
            end
        end
        setdirty!(state_orig)
        state_prin = deepcopy(state_orig)
        principal_value!(state_prin)
        for joint_k = tree_joints(state_orig.mechanism)
            joint_type_k = joint_type(joint_k)
            q_orig = state_orig.q[joint_k]
            q_prin = state_prin.q[joint_k]
            if joint_type_k isa SPQuatFloating
                rot_orig = rotation(joint_type_k, q_orig)
                rot_prin = rotation(joint_type_k, q_prin)
                @test isapprox(rot_orig, rot_prin)
                @test (rot_prin.x^2 + rot_prin.y^2 + rot_prin.z^2) < (1.0 + 1.0e-14)
                @test (1.0 + 1.0e-14) < (rot_orig.x^2 + rot_orig.y^2 + rot_orig.z^2)
            elseif joint_type_k isa QuaternionFloating || joint_type_k isa QuaternionSpherical
                rot_orig = rotation(joint_type_k, q_orig)
                rot_prin = rotation(joint_type_k, q_prin)
                @test isapprox(rot_orig, rot_prin)
                @test rot_prin.w > 0
            else
                @test q_orig == q_prin
            end
            if joint_type_k isa QuaternionFloating
                @test translation(joint_type_k, q_orig) === translation(joint_type_k, q_prin)
            end
        end
    end
end


text:
using Plots
using Random
using Pkg; Pkg.add("LaTeXStrings")
using LaTeXStrings

# Set seed
Random.seed!(123)

# Generate MC points
n = 2000
points_x = rand(n)
points_y = rand(n)

# Label points in circle
labels = [x^2 + y^2 <= 1 for (x,y) in zip(points_x, points_y)]

# Estaimte of pi
π_estimate = 4*sum(labels)/n

# Create plot of MC points
plot(points_x, points_y, seriestype = :scatter,
     group = .!labels,
     title = L"\textrm{Monte Carlo } \pi \textrm{ estimation}",
     legend = false,
     bg = RGB(247/255, 236/255, 226/255))

# Add semi-circle quadrant plot
circle_x = 0:0.01:1
circle_y = sqrt.(1 .- circle_x.^2)
plot!(circle_x, circle_y, seriestype= :line, legend = false, lc = :yellow)

# Generate plots comparing π estimate across multiple n
n_values = [5000, 10000, 100000, 1e6, 2e6]
π_estimates_error = zeros(5)

for i in 1:5
    n_value = Int64(n_values[i])
    points_x = rand(n_value)
    points_y = rand(n_value)
    estimate = 4*sum([x^2 + y^2 <= 1 for (x,y) in zip(points_x, points_y)])/
        n_value
    π_estimates_error[i] = 100*abs(π - estimate)/π
end

println(π_estimates_error)

# savefig("../images/pi-monte-carlo.svg")

# Plot area under circle curve
plot(circle_x, circle_y, legend = false,
    fillrange = [zeros(length(circle_x)), circle_y],
    alpha = .5,
    lc = :yellow)

ext: jl

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
# Copyright 2021 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Executes the identifiability experiement with some instances of the netlib
# database.
# To run this script execute
# $ julia --project run_netlib_infeasibility_experiment.jl <path_where_to_store_results> <path_to_directory_containing_mps_files>
#
# By default the script runs a specific set of twelve instances. The script also
# optionally accepts a specific instance to run, for example,
# $ julia --project run_netlib_infeasibility_experiment.jl <path_where_to_store_results> <path_to_directory_containing_mps_files> box1

import FirstOrderLp

include("primal_dual_hybrid_gradient.jl")
include("utils.jl")

function run_netlib_instance(instance_name::String)
  params = PdhgParameters(
    nothing, # step_size (forces the solver to use a provably correct step size)
    1.0, # primal_weight
    1000, # printing frequency
    false, # verbose
    1000000, # iteration limit
    nothing, # two_dimensional_subspace
    1e-8, # eps_primal_infeasible
    1e-8, # eps_dual_infeasible
    false, # store_iterates
    nothing, # reference_point
  )
  print("About to start solving " * instance_name * "\n")
  instance_path = joinpath(ARGS[2], (instance_name * ".mps"))
  lp = FirstOrderLp.qps_reader_to_standard_form(instance_path)
  output = optimize(params, lp)
  generate_infeasibility_error_plot(
    output,
    joinpath(ARGS[1], instance_name * "-error.pdf"),
  )
end

function main()
  if length(ARGS) == 3
    run_netlib_instance(ARGS[3])
  else
    for instance_name in [
      "box1",
      "galenet",
      "itest2",
      "itest6",
      "woodinfe",
      "ex72a",
      "ex73a",
      "bgdbg1",
      "bgindy",
      "cplex1",
      "chemcom",
      "mondou2",
    ]
      run_netlib_instance(instance_name)
    end
  end
end

main()


text:
# test correlation

T = Float32
fs = 100.  # sampling rate
N = 2^14   # number of points
Nwin = 10
cc_len = 40.96 # length of noise window
cc_step = 40.96 # step between windows
freqmin = 10. # minimum frequency
freqmax = 20. # maximum frequency
corners = 4   # number of corners
zerophase = true # do zerophase filtering
maxlag = 20. # seconds of coda to save
starttime = d2u(DateTime(Date(now()))) # starttime for SeisChanel / SeisData
rand_net()    = uppercase(String(rand('A':'Z', 2)))
rand_sta()    = String(rand('A':'Z', rand(3:5)))
rand_loc()    = rand() < 0.5 ? "" : lpad(rand(0:9), 2, "0")
rand_id() = string(rand_net(), ".", rand_sta(), ".", rand_loc(), ".HH", rand('A':'Z'))

@testset "clean up" begin
    A = rand(T,N,Nwin)
    B = clean_up(A,freqmin,freqmax,fs)
    C = clean_up(A,Int(freqmin),Int(freqmax),Int(fs))
    @test size(B) == size(A)
    @test eltype(B) == eltype(A)
    @test eltype(C) == eltype(A)
    @test B != A
    @test C == B

    # test in-place
    clean_up!(A,freqmin,freqmax,fs)
    @test A == B
    @test A == C

    # test RawData
    R = RawData()
    R.x = rand(T,N,Nwin)
    R.fs = fs
    R.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    Rclean = clean_up(R,freqmin,freqmax)
    RInt = clean_up(R,Int(freqmin),Int(freqmax))
    @test size(Rclean.x) == size(R.x)
    @test eltype(Rclean.x) == eltype(R.x)
    @test size(RInt.x) == size(R.x)
    @test eltype(RInt.x) == eltype(R.x)
    @test Rclean.x != R.x
    @test Rclean.freqmin == freqmin
    @test Rclean.freqmax == freqmax
    @test RInt.freqmin == freqmin
    @test RInt.freqmax == freqmax

    # test in-place
    clean_up!(R,freqmin,freqmax)
    @test R.x == Rclean.x
    @test R.x == RInt.x
    @test R.freqmin == freqmin
    @test R.freqmax == freqmax

    # test CorrData
    C = CorrData()
    C.corr = rand(T,N,Nwin)
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    Cclean = clean_up(C,freqmin,freqmax)
    CInt = clean_up(C,Int(freqmin),Int(freqmax))
    @test size(Cclean.corr) == size(C.corr)
    @test eltype(Cclean.corr) == eltype(C.corr)
    @test Cclean.corr != C.corr
    @test Cclean.freqmin == freqmin
    @test Cclean.freqmax == freqmax
    @test CInt.freqmin == freqmin
    @test CInt.freqmax == freqmax

    # test in-place
    clean_up!(C,freqmin,freqmax)
    @test C.corr == Cclean.corr
    @test C.corr == CInt.corr
    @test C.freqmin == freqmin
    @test C.freqmax == freqmax
end

@testset "correlation" begin
    # test correlate function on arrays
    A = rand(T,N,Nwin) .- T(0.5)
    FFT = rfft(A,1)

    C = correlate(FFT,FFT,N,Int(maxlag * fs))
    @test size(C) == (Int(2 * maxlag * fs) + 1,Nwin)
    @test eltype(C) == eltype(A)

    # test autocorrelation peak occurs at zero lag
    t = -maxlag:1/fs:maxlag
    maxinds = argmax(C,dims=1)
    @test all([t[maxinds[ii][1]] == 0 for ii in length(maxinds)])

    # test phase correlation on array
    A = rand(T,N,Nwin) .- T(0.5)
    FFT = fft(A,1)

    C = phasecorrelate(FFT,FFT,N,Int(maxlag * fs))
    @test size(C) == (Int(2 * maxlag * fs) + 1,Nwin)
    @test eltype(C) == eltype(A)

    # test autocorrelation peak occurs at zero lag
    t = -maxlag:1/fs:maxlag
    maxinds = argmax(C,dims=1)
    @test all([t[maxinds[ii][1]] == 0 for ii in length(maxinds)])

    # test correlation with FFTData
    Ch = SeisIO.RandSeis.randSeisChannel(c=false,s=true)
    ungap!(Ch)
    Ch.x = rand(T,Int(cc_len*Nwin*fs +1)) .- T(0.5)
    Ch.fs = fs
    Ch.t = [1 starttime * 1e6;length(Ch.x) 0]
    Ch.loc = GeoLoc()
    R = RawData(Ch,cc_len,cc_step)
    F = rfft(R)
    C = correlate(F,F,maxlag)
    CInt = correlate(F,F,Int(maxlag))
    @test isa(C,CorrData) # test return type
    @test C.comp == repeat(Ch.id[end],2) # test component name
    @test size(C.corr) == (Int(2 * maxlag * fs) + 1,Nwin) # test size
    @test eltype(C.corr) == eltype(Ch.x) # test type stability
    maxinds = argmax(C.corr,dims=1)
    @test all([t[maxinds[ii][1]] == 0 for ii in length(maxinds)]) # test max args
    @test C.rotated == false # test rotation
    @test C.corr_type == "CC" # test cross-correlation
    @test C.maxlag == CInt.maxlag
    @test C.corr == CInt.corr

    # test windows that do not overlap
    F1 = deepcopy(F)
    F2 = deepcopy(F)
    F2.t[end] = rand(eltype(F2.t))
    C1 = correlate(F1,F2,maxlag)
    @test size(C1.corr) == (Int(2 * maxlag * fs) + 1,Nwin-1) # test size

    # test with no matching windows
    F2.t = rand(eltype(F2.t),size(F2.fft,2))
    @test_throws ArgumentError correlate(F1,F2,maxlag)

    # test phase correlation
    P = phase(R)
    PCC = correlate(P,P,maxlag,corr_type="PCC")
    @test isa(PCC,CorrData) # test return type
    @test PCC.comp == repeat(Ch.id[end],2) # test component name
    @test size(PCC.corr) == (Int(2 * maxlag * fs) + 1,Nwin) # test size
    @test eltype(PCC.corr) == eltype(Ch.x) # test type stability
    maxinds = argmax(PCC.corr,dims=1)
    @test all([t[maxinds[ii][1]] == 0 for ii in length(maxinds)]) # test max args
    @test PCC.rotated == false # test rotation
    @test PCC.corr_type == "PCC" # test cross-correlation

    # test unrecognized type of correlation
    @test_throws ArgumentError correlate(F,F,maxlag,corr_type="AUTO")

end

@testset "whitening" begin

    pad = 50
    padfreq = pad / fs
    # test whitening array
    A = rand(T,N,Nwin) .- T(0.5)
    F = rfft(A,1)
    Fwhite = whiten(F,freqmin,freqmax,fs,N,pad=pad)

    # test whitening
    freqs = FFTW.rfftfreq(N,fs)
    lowband = findall(0 + eps(T) .< freqs .< freqmin - padfreq)
    passband = findall(freqmin .<= freqs .<= freqmax)
    highband = findall(freqs .> freqmax + padfreq)

    # test filtering
    @test all(Fwhite[lowband,:] .== 0)
    @test all(Fwhite[passband,:] .!= 0)
    @test all(Fwhite[highband,:] .== 0)
    @test size(Fwhite) == size(F)
    @test eltype(Fwhite) == eltype(F)

    # test with low whitening passband
    Flow = whiten(F,0.,freqmax,fs,N,pad=pad)
    @test all(abs.(Flow[1:pad,:]) .< 1.)

    # test with high whitening passband
    Fhigh = whiten(F,freqmin,fs/2,fs,N,pad=pad)
    @test all(abs.(Fhigh[end-pad:end,:]) .< 1.)

    # test with integer arguments
    FInt = whiten(F,Int(freqmin),Int(freqmax),Int(fs),N,pad=pad)

    # test in-place
    whiten!(F,freqmin,freqmax,fs,N,pad=pad)
    @test F == Fwhite
    @test F == FInt

    # test FFTData
    Ch = SeisIO.RandSeis.randSeisChannel(c=false,s=true)
    Ch.loc = GeoLoc()
    ungap!(Ch)
    Ch.x = rand(T,Int(cc_len*Nwin*fs +1)) .- T(0.5)
    Ch.fs = fs
    Ch.t = [1 starttime * 1e6;length(Ch.x) 0]
    R = RawData(Ch,cc_len,cc_step)
    F = rfft(R)

    Fwhite = whiten(F,freqmin,freqmax)
    @test isa(Fwhite,FFTData)
    @test size(Fwhite.fft) == size(F.fft)
    @test Fwhite.whitened == true
    @test Fwhite.freqmin == freqmin
    @test Fwhite.freqmax == freqmax
    @test Fwhite.fft != F.fft

    # test low/high frequency
    Flowhigh = whiten(F,0.,fs)
    @test Flowhigh.whitened == true
    @test Flowhigh.freqmin == F.freqmin
    @test Flowhigh.freqmax == F.freqmax

    # test low frequency
    Flow = whiten(F,0.,freqmax)
    @test Flow.whitened == true
    @test Flow.freqmin == F.freqmin
    @test Flow.freqmax == freqmax

    # test low frequency
    Fhigh = whiten(F,freqmin,fs)
    @test Fhigh.whitened == true
    @test Fhigh.freqmin == freqmin
    @test Fhigh.freqmax == F.freqmax

    # test with integer input
    FInt = whiten(F,Int(freqmin),Int(freqmax))
    @test FInt.freqmin == freqmin
    @test FInt.freqmax == freqmax

    # test in-place
    whiten!(F,freqmin,freqmax)
    @test F.fft == Fwhite.fft
    @test F.whitened == true
    @test F.freqmin == freqmin
    @test F.freqmax == freqmax
    @test FInt == F

    # test RawData
    Rwhite = whiten(R,freqmin,freqmax)
    @test isa(Rwhite,RawData)
    @test size(Rwhite.x) == size(R.x)
    @test Rwhite.whitened == true
    @test Rwhite.freqmin == freqmin
    @test Rwhite.freqmax == freqmax
    @test Rwhite.x != R.x

    # test low/high frequency
    Rlowhigh = whiten(R,0.,fs)
    @test Rlowhigh.whitened == true
    @test Rlowhigh.freqmin == R.freqmin
    @test Rlowhigh.freqmax == R.freqmax

    # test low frequency
    Rlow = whiten(R,0.,freqmax)
    @test Rlow.whitened == true
    @test Rlow.freqmin == R.freqmin
    @test Rlow.freqmax == freqmax

    # test low frequency
    Rhigh = whiten(R,freqmin,fs)
    @test Rhigh.whitened == true
    @test Rhigh.freqmin == freqmin
    @test Rhigh.freqmax == R.freqmax

    # test with integer input
    RInt = whiten(R,Int(freqmin),Int(freqmax))
    @test RInt.freqmin == freqmin
    @test RInt.freqmax == freqmax

    # test in-place
    whiten!(R,freqmin,freqmax)
    @test R.x == Rwhite.x
    @test R.whitened == true
    @test R.freqmin == freqmin
    @test R.freqmax == freqmax
    @test RInt == R

    # test coherence with smoothing
    R = RawData(Ch,cc_len,cc_step)
    F = rfft(R)
    Fcoh = coherence(F,20,0.01)
    @test isa(Fcoh,FFTData)
    @test size(Fcoh.fft) == size(F.fft)
    @test eltype(Fcoh.fft) == eltype(F.fft)
    @test Fcoh.fft != F.fft

    # test coherence no smoothing
    Fcoh = coherence(F,20)
    @test isa(Fcoh,FFTData)
    @test size(Fcoh.fft) == size(F.fft)
    @test eltype(Fcoh.fft) == eltype(F.fft)
    @test Fcoh.fft != F.fft

    # test inplace
    coherence!(F,20)
    @test F.fft == Fcoh.fft

    # test deconvolution with smoothing
    R = RawData(Ch,cc_len,cc_step)
    F = rfft(R)
    Fdec = deconvolution(F,20,0.01)
    @test isa(Fdec,FFTData)
    @test size(Fdec.fft) == size(F.fft)
    @test eltype(Fdec.fft) == eltype(F.fft)
    @test Fdec.fft != F.fft

    # test deconvolution no smoothing
    Fdec = deconvolution(F,20)
    @test isa(Fdec,FFTData)
    @test size(Fdec.fft) == size(F.fft)
    @test eltype(Fdec.fft) == eltype(F.fft)
    @test Fdec.fft != F.fft

    # test inplace
    deconvolution!(F,20)
    @test F.fft == Fdec.fft
end


text:
using DiffEqSensitivity, QuasiMonteCarlo, Test

function ishi_batch(X)
    A= 7
    B= 0.1
    @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
end
function ishi(X)
    A= 7
    B= 0.1
    sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
end

function linear_batch(X)
    A= 7
    B= 0.1
    @. A*X[1,:]+B*X[2,:]
end
function linear(X)
    A= 7
    B= 0.1
    A*X[1]+B*X[2]
end

n = 600000
lb = -ones(4)*π
ub = ones(4)*π
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(n,lb,ub,sampler)

res1 = gsa(ishi,Sobol(),A,B)
res2 = gsa(ishi_batch,Sobol(),A,B,batch=true)

@test res1.S1 ≈ [0.3139335358797363, 0.44235918402206326, 0.0, 0.0] atol=1e-4
@test res2.S1 ≈ [0.3139335358797363, 0.44235918402206326, 0.0, 0.0] atol=1e-4

@test res1.ST ≈ [0.5576009081644232, 0.44237102330046346, 0.24366241588532553, 0.0] atol=1e-4
@test res2.ST ≈ [0.5576009081644232, 0.44237102330046346, 0.24366241588532553, 0.0] atol=1e-4

res1 = gsa(linear,Sobol(),A,B)
res2 = gsa(linear_batch,Sobol(),A,B,batch=true)

@test res1.S1 ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4
@test res2.S1 ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4

@test res1.ST ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4
@test res2.ST ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4

#=
library(sensitivity)
ishigami.fun <- function(X) {
  A <- 7
  B <- 0.1
  sin(X[, 1]) + A * sin(X[, 2])^2 + B * X[, 3]^4 * sin(X[, 1])
}
n <- 600000
X1 <- data.frame(matrix(runif(4 * n,-pi,pi), nrow = n))
X2 <- data.frame(matrix(runif(4 * n,-pi,pi), nrow = n))
sobol2007(ishigami.fun, X1, X2)
sobolSalt(ishigami.fun, X1, X2, scheme="A")
sobolSalt(ishigami.fun, X1, X2, scheme="B")

library(sensitivity)
ishigami.fun <- function(X) {
  A <- 7
  B <- 0.1
  A * X[, 1] + B * X[, 2]
}
n <- 6000000
X1 <- data.frame(matrix(runif(4 * n,-pi,pi), nrow = n))
X2 <- data.frame(matrix(runif(4 * n,-pi,pi), nrow = n))
sobol2007(ishigami.fun, X1, X2)
sobolSalt(ishigami.fun, X1, X2, scheme="A")
sobolSalt(ishigami.fun, X1, X2, scheme="B")
=#

function ishi_linear(X)
    A= 7
    B= 0.1
    [sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1]),A*X[1]+B*X[2]]
end

function ishi_linear_batch(X)
    A= 7
    B= 0.1
    X1 = @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
    X2 = @. A*X[1,:]+B*X[2,:]
    vcat(X1',X2')
end

res1 = gsa(ishi_linear,Sobol(),A,B)
res2 = gsa(ishi_linear_batch,Sobol(),A,B,batch=true)

# Now both tests together

@test res1.S1[1,:] ≈ [0.3139335358797363, 0.44235918402206326, 0.0, 0.0] atol=1e-4
@test res2.S1[1,:] ≈ [0.3139335358797363, 0.44235918402206326, 0.0, 0.0] atol=1e-4

@test res1.ST[1,:] ≈ [0.5576009081644232, 0.44237102330046346, 0.24366241588532553, 0.0] atol=1e-4
@test res2.ST[1,:] ≈ [0.5576009081644232, 0.44237102330046346, 0.24366241588532553, 0.0] atol=1e-4

@test res1.S1[2,:] ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4
@test res2.S1[2,:] ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4

@test res1.ST[2,:] ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4
@test res2.ST[2,:] ≈ [0.9997953478183109, 0.0002040399766938839, 0.0, 0.0] atol=1e-4


text:
using Dates: value

function Base.print(io::IO, tz::FixedTimeZone)
    if get(io, :compact, true)
        isempty(tz.name) ? print(io, "UTC", tz.offset) : print(io, tz.name)
    else
        offset_str = "UTC" * offset_string(tz.offset, true)  # Use ISO 8601 for comparision
        if isempty(tz.name)
            print(io, offset_str)
        elseif tz.name != offset_str && !(value(tz.offset) == 0 && tz.name in ("UTC", "GMT"))
            print(io, tz.name, " (UTC", offset_string(tz.offset), ")")
        else
            print(io, tz.name)
        end
    end
end

function Base.print(io::IO, tz::VariableTimeZone)
    if get(io, :compact, true)
        print(io, tz.name)
    else
        trans = tz.transitions

        # Retrieve the "modern" time zone transitions. We'll treat the latest transitions as
        # the same as the transitions for `now()` since these future transitions should be
        # based upon the same rules.
        if tz.cutoff === nothing || length(trans) == 1
            trans = trans[end:end]
        else
            trans = trans[end-1:end]

            # Attempt to show a standard time offset before daylight saving time offset.
            # Sorting should work as long as the DST adjustment is always positive. Fixes
            # differences between the north and south hemispheres.
            sort!(trans, by=el -> el.zone.offset)
        end

        # Show standard time offset before daylight saving time offset.
        print(
            io,
            tz.name,
            " (", join(["UTC" * offset_string(t.zone.offset) for t in trans], "/"), ")",
        )
    end
end

function Base.print(io::IO, t::Transition)
    print(io, t.utc_datetime, " ")
    show(io, MIME("text/plain"), t.zone.offset)  # Long-form
    !isempty(t.zone.name) && print(io, " (", t.zone.name, ")")
end

Base.print(io::IO, zdt::ZonedDateTime) = print(io, DateTime(zdt, Local), zdt.zone.offset)

function Base.show(io::IO, tz::FixedTimeZone)
    if istimezone(tz.name, Class(:ALL)) && isequal(tz, TimeZone(tz.name, Class(:ALL)))
        print(io, "tz\"$(tz.name)\"")
    else
        std = Dates.value(tz.offset.std)
        dst = Dates.value(tz.offset.dst)

        params = [repr(tz.name), repr(std)]
        dst != 0 && push!(params, repr(dst))
        print(io, FixedTimeZone, "(", join(params, ", "), ")")
    end
end

function Base.show(io::IO, tz::VariableTimeZone)
    # Compat printing when the time zone can be constructed with `@tz_str`
    if istimezone(tz.name, Class(:ALL)) && isequal(tz, TimeZone(tz.name, Class(:ALL)))
        print(io, "tz\"$(tz.name)\"")

    # Compact printing of a custom time zone which is non-constructable
    elseif get(io, :compact, false)
        print(io, VariableTimeZone, "(")
        show(io, tz.name)
        print(io, ", ...)")

    # Verbose printing which should print a fully constructable `VariableTimeZone`.
    else
        # Force `:compact => false` to make the force the transition vector printing into
        # long form.
        print(io, VariableTimeZone, "(")
        show(io, tz.name)
        print(io, ", ")
        show(IOContext(io, :compact => false), tz.transitions)
        print(io, ", ")
        show(io, tz.cutoff)
        print(io, ")")
    end
end

function Base.show(io::IO, t::Transition)
    if get(io, :compact, false) || get(io, :compact_el, false)
        print(io, t)
    else
        # Fallback to calling the default show instead of reimplementing it.
        invoke(show, Tuple{IO, Any}, io, t)
    end
end

function Base.show(io::IO, zdt::ZonedDateTime)
    if VERSION <= v"1.5.0-DEV.163" && get(io, :compact_el, false)
        print(io, zdt)
    else
        values = [
            yearmonthday(zdt)...
            hour(zdt)
            minute(zdt)
            second(zdt)
            millisecond(zdt)
        ]
        index = something(findlast(!iszero, values), 1)
        params = [
            map(repr, values[1:index]);
            repr(timezone(zdt); context=:compact => true)
        ]

        print(io, ZonedDateTime, "(", join(params, ", "), ")")
    end
end

Base.show(io::IO, ::MIME"text/plain", t::Transition) = print(io, t)
Base.show(io::IO, ::MIME"text/plain", tz::TimeZone) = print(IOContext(io, :compact => false), tz)
Base.show(io::IO, ::MIME"text/plain", zdt::ZonedDateTime) = print(io, zdt)

# https://github.com/JuliaLang/julia/pull/33290
if VERSION >= v"1.5.0-DEV.224"
    Base.typeinfo_implicit(::Type{ZonedDateTime}) = true
end

# Use compact printing on certain element types for Julia versions before:
# https://github.com/JuliaLang/julia/pull/34387
if VERSION <= v"1.5.0-DEV.163"
    function Base.show(io::IO, m::MIME"text/plain", X::AbstractArray{T}) where T <: Union{Transition, ZonedDateTime}
        # Specify a custom IO context which we can use to perform compact printing of
        # elements.
        io = IOContext(io, :compact_el => true)
        invoke(show, Tuple{IO, MIME"text/plain", AbstractArray}, io, m, X)
    end
end


text:
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
using ..Units
using Parameters, Unitful

# Computational routines needed for computations with the MaterialParams structure

# with tuple & vector - apply for all phases in MatParam
function compute_param!(fn::F, rho::Vector{_T}, MatParam::NTuple{N,AbstractMaterialParamsStruct}, P::_T=zero(_T),T::_T=zero(_T)) where {F,N,_T}
    rho .= map(x->fn(x,P,T), MatParam)
end

# each individual calcuation
function compute_param(fn::F, MatParam::NTuple{N,AbstractMaterialParamsStruct}, P::_T=zero(_T),T::_T=zero(_T)) where {F,N,_T}
    map(x->fn(x,P,T), MatParam)
end

#---------------------------------------------------------------------------------------------------------------------------#
#Computational routines for Phases

#performs computation given a single Phase
function compute_param(fn::F, MatParam::NTuple{N,AbstractMaterialParamsStruct}, Phase::Int64, P::_T=zero(_T),T::_T=zero(_T)) where {F,N,_T}
    Phase_tup = ntuple(i->MatParam[i].Phase, Val(N))
    ind = find_ind(Phase_tup, Phase)
    return fn(MatParam, P, T)[ind]
end

function compute_param!(fn::F, rho::AbstractArray{_T, ndim}, MatParam::NTuple{N,AbstractMaterialParamsStruct}, Phases::AbstractArray{_I, ndim}, P=nothing, T=nothing) where {F,ndim,N,_T,_I<:Integer}
    Phase_tup = ntuple(i->MatParam[i].Phase, Val(N))

    # check if we compute with P/T or use default values
    isnothing(T) ? compute_T = false : compute_T = true
    isnothing(P) ? compute_P = false : compute_P = true

    Tval = zero(_T)
    Pval = zero(_T)

    @inbounds for I in eachindex(Phases)
        phase   = find_ind(Phase_tup, Phases[I])

        # Extract relevant value if requested
        if compute_T; Tval = T[I]; end
        if compute_P; Pval = P[I]; end

        # this computes density for ALL phases, which is a bit of an overkill as we only need density for a single phase
        rho_tup = fn(MatParam, Pval, Tval)
        rho[I]  = rho_tup[phase]
    end
end

#Phase PhaseRatios

function compute_param!(fn::F, rho::AbstractArray{_T, N}, MatParam::NTuple{K,AbstractMaterialParamsStruct}, PhaseRatios::AbstractArray{_T, M}, P=nothing, T=nothing) where {F,_T<:AbstractFloat, N,M, K}
    if M!=(N+1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    # check if we compute with P/T or use default values
    isnothing(T) ? compute_T = false : compute_T = true
    isnothing(P) ? compute_P = false : compute_P = true

    Tval = zero(_T)
    Pval = zero(_T)
    if compute_P
        Rindex = CartesianIndices(P)
    else
        Rindex = CartesianIndices(T)
    end

    @inbounds for I in Rindex
        frac    = view(PhaseRatios, Tuple(I)..., 1:K)    # fraction of each phase @ point I

        # Extract relevant value if requested
        if compute_T; Tval = T[I]; end
        if compute_P; Pval = P[I]; end

        # compute point-wise density:
        rho[I] = compute_param_times_frac(fn, frac, MatParam, Pval, Tval)
    end

    return
end

#Multiplies parameter with the fraction of a phase
function compute_param_times_frac(fn::F, PhaseRatios::AbstractArray{_T, 1}, MatParam::NTuple{N,AbstractMaterialParamsStruct}, P::_T, T::_T) where {F,_T, N}

    value_tup = fn(MatParam, P, T)

    # Sum & multiply density with fraction
    val = zero(_T)
    @inbounds for j=1:N
        val += PhaseRatios[j]*value_tup[j]
    end

    return val
end


text:
using ProgressMeter
using .Threads
using Lasso
using NPZ
using Statistics
using Random
using StatsBase
using LinearAlgebra

struct InvResult{T,S}
    i::Int
    λmin::T
    dn::T
    m::S
end

function fit_model(G,dv,i)
    dn = std(dv)
    dv ./= dn
    path = fit(LassoPath, G, dv, λminratio = 1e-5, standardize=false, intercept=false, maxncoef=size(G,2))
    λmin = path.λ[argmin(aicc(path))]
    m = coef(path; select=MinAICc())
    return InvResult(i, λmin, dn, m)
end

T = Float64

Gzigzag = npzread("Kernels/G_zigzag.npy")
Gspiral = npzread("Kernels/G_spiral.npy")
Gcross = npzread("Kernels/G_cross.npy")
Grandom = npzread("Kernels/G_random.npy")

dzigzag = npzread("Testing/zigzag.npy")
dspiral = npzread("Testing/spiral.npy")
dcrossing = npzread("Testing/crossing.npy")
drandom = npzread("Testing/random.npy")

p = Progress(4*size(dzigzag,2), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

for (G, d, name) in zip([Gzigzag, Gspiral, Gcross, Grandom], [dzigzag, dspiral, dcrossing, drandom], ["zigzag", "spiral", "crossing", "random"])
        results = Array{Any,1}(undef, size(d,2))
        uresults = Array{Any,1}(undef, size(d,2))

        G = G / std(G)

        @threads for j = 1:size(d,2)
            results[j] = fit_model(G, d[:,j], j)
            next!(p)
        end

        sres = reduce(hcat, [x.m.*x.dn for x in results])

        npzwrite("Testing/$(name)_res.npy", Array{Float64}(sres))

end

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
###############################################################################
#### Single spin ##############################################################
###############################################################################

# SingleSpin(x::T, y::T, z::T) where T = [x, y, z]

###############################################################################
#### Product states ###########################################################
###############################################################################

ProductState(states::Array{Array{Float64, 1}, 1}) where T = hcat(states...)
ProductState(N::Int, state::AbstractVector{Float64}) where T = ProductState([state for i in 1:N])

###############################################################################
#### magnetization ############################################################
###############################################################################

function magnetization(s::AbstractMatrix, x::Symbol=:x)
    match = Dict(:x => 1, :y => 2, :z=>3)

    N = size(s, 2)
    1/N * sum(s[match[x], :])
end

magnetization(sol::OrdinaryDiffEq.ODECompositeSolution, t) = magnetization(sol(t)')
magnetization(sim::EnsembleSolution, t) = magnetization(timepoint_mean(sim, t)')
###############################################################################
#### dtwa states ############################################################
###############################################################################

function dtwa_state(x::AbstractVector)
    res = normalize(x)
    nullspace(permutedims(x))[:, 1]
    vec(res .+ sum(rand([-1, 1], 1, 2)  .* nullspace(permutedims(res)), dims=2))
end

function dtwa_state(x::AbstractMatrix)
    res = similar(x)
    for i in 1:size(x, 2)
        res[:, i] .= dtwa_state(x[:, i])
    end
    res
end

# end #module SpinState

max_stars_repo_name: titusfranz/DTWA

max_issues_repo_name: titusfranz/DTWA

max_forks_repo_name: titusfranz/DTWA


text:
module RCAM
#Residual Constrained Alternating Minimization

    using Distributed
    using DistributedArrays
    using LinearAlgebra
    include("common.jl")
    include("dclr.jl")
    include("dclrMN.jl")

end


text:
################################################################################
# Note: The following is based on the nomenclature in the paper
# "Revised Release on the IAPWS Formulation 1995 for the Thermodynamic
# Properties of Ordinary Water Substance for General and Scientific Use"
# by IAPWS, the International Association for the Properties of Water and Steam.
#
# It is more or less as follows:
#
# ϕ : Dimensionless Helmholtz energy of water as a fucntion of density and temperature
# δ : Dimensionless density of Water
# τ : Dimensionless temperature of water
# n, d, t, α, β, γ, ϵ, A, B, C, D : these are parameters used in fitting the
# equation of state to data. As far as we're concerned, they are constants with
# no real physical significance.
#
# Note that although the code was written by consulting the IAPWS paper, it is
# applicable to any
################################################################################

struct IAPWSTypeMultiparameterEOS{T_} <: AbstractEquationOfState where T_ <: AbstractFloat
    Tc::T_
    ρc::T_
    pc::T_
    ω::T_
    R::T_

    #parameters for ideal gas portion
    n₀::Vector{T_}
    γ₀::Vector{T_}

    # parameters for residual term 1
    n₁::Vector{T_}
    d₁::Vector{T_}
    t₁::Vector{T_}

    # parameters for residual term 2
    n₂::Vector{T_}
    d₂::Vector{T_}
    t₂::Vector{T_}
    c₂::Vector{T_}

    # parameters for residual term 3
    n₃::Vector{T_}
    d₃::Vector{T_}
    t₃::Vector{T_}
    α₃::Vector{T_}
    β₃::Vector{T_}
    γ₃::Vector{T_}
    ϵ₃::Vector{T_}

    # parameters for residual term 4
    n₄::Vector{T_}
    a₄::Vector{T_}
    b₄::Vector{T_}
    β₄::Vector{T_}
    A₄::Vector{T_}
    B₄::Vector{T_}
    C₄::Vector{T_}
    D₄::Vector{T_}

end

function ϕ₀(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    first = log(δ) + eos.n₀[1] + eos.n₀[2]*τ + eos.n₀[3]*log(τ)
    end_ = sum(eos.n₀[4:8].*log.(1-exp.(-eos.γ₀.*τ)))
    return first + end_
end

function res_1(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum(eos.n₁ .* (δ .^ eos.d₁) .* (τ .^ eos.t₁))
end

function res_2(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum(eos.n₂ .* (δ .^ eos.d₂) .* (τ .^ eos.t₂) .* exp.(-δ .^ eos.c₂))
end

function res_3(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum(eos.n₃ .* (δ .^ eos.d₃) .* (τ .^ eos.t₃) .* exp.(-eos.α₃ .* (δ - eos.ϵ₃).^2 - eos.β₃ .* (τ - eos.γ₃).^2))
end

function θ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum((1-τ) + eos.A₄.*((δ - 1).^2).^(1 ./(2 .* eos.β₄)))
end

function Δ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum(θ(eos, δ, τ).^2 + eos.B₄ .* ((δ - 1).^2) .* eos.a₄)
end

function Ψ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum(exp(-eos.C₄ .* (δ - 1)^2 - eos.D₄ .* (τ - 1)^2))
end

function res_4(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return sum(eos.n₄ .* Δ(eos,δ,τ).^eos.b₄ .* δ * Ψ(eos,δ,τ))
end

function ϕr(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return res_1(eos, δ, τ)+res_2(eos, δ, τ)+res_3(eos, δ, τ)+res_4(eos, δ, τ)
end

function ϕ₀_δ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(δ_) = ϕ₀(eos, δ_, τ)
    return ForwardDiff.derivative(f, δ)
end
function ϕ₀_δδ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(δ_) = ϕ₀_δ(eos, δ_, τ)
    return ForwardDiff.derivative(f, δ)
end
function ϕ₀_τ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(τ_) = ϕ₀(eos, δ, τ_)
    return ForwardDiff.derivative(f, τ)
end
function ϕ₀_ττ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(τ_) = ϕ₀_τ(eos, δ, τ_)
    return ForwardDiff.derivative(f, τ)
end
function ϕ₀_δτ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    return 0.0
end

function ϕr_δ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(δ_) = ϕr(eos, δ_, τ)
    return ForwardDiff.derivative(f, δ)
end
function ϕr_δδ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(δ_) = ϕr_δ(eos, δ_, τ)
    return ForwardDiff.derivative(f, δ)
end
function ϕr_τ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(τ_) = ϕr(eos, δ, τ_)
    return ForwardDiff.derivative(f, τ)
end
function ϕr_ττ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(τ_) = ϕr_τ(eos, δ, τ_)
    return ForwardDiff.derivative(f, τ)
end
function ϕr_δτ(eos::IAPWSTypeMultiparameterEOS, δ, τ)
    f(τ_) = ϕr_δ(eos, δ, τ_)
    return ForwardDiff.derivative(f, τ)
end

function p(eos::IAPWSTypeMultiparameterEOS, ρ, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    #@printf "δ = %f\n" δ
    #@printf "τ = %f\n" τ
    #@printf "derivative is %f\n" ϕr_δ(eos, δ, τ)
    p_nondim = 1 + δ*ϕr_δ(eos, δ, τ)
    return p_nondim * ρ*eos.R*T
end

function ρ_from_p(eos::IAPWSTypeMultiparameterEOS, p_, T)
    function solveit(ρ_guess)
        return abs(p_ - p(eos, ρ_guess, T))
    end
    a = optimize(solveit, 0.00001, 1500.0)
    print(a)
end

function u(eos::IAPWSTypeMultiparameterEOS, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    u_nondim = τ * (ϕ₀_τ(eos, δ, τ) + ϕr_τ(eos, δ, τ))
    return u_nondim * T*eos.R
end

function s(eos::IAPWSTypeMultiparameterEOS, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    s_nondim = τ * (ϕ₀_τ(eos, δ, τ) + ϕr_τ(eos, δ, τ)) - (ϕ₀(eos, δ, τ) + ϕr(eos, δ, τ))
    return s_nondim * eos.R
end

function H(eos::IAPWSTypeMultiparameterEOS, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    H_nondim = 1 + τ * (ϕ₀_τ(eos, δ, τ) + ϕr_τ(eos, δ, τ)) + δ*ϕr_δ(eos, δ, τ)
    return H_nondim * T*eos.R
end

function Cp(eos::IAPWSTypeMultiparameterEOS, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    part1 = -τ^2 * (ϕ₀_ττ(eos, δ, τ) + ϕr_ττ(eos, δ, τ))
    p2_top = (1 + δ*ϕr_δ(eos, δ, τ) - δ*τ*ϕr_δτ(eos, δ, τ))^2
    p2_bot = 1 + 2δ*ϕr_δ(eos, δ, τ) + δ^2*ϕr_δδ(eos, δ, τ)
    Cp_nondim = part1 + p2_top/p2_bot
    return Cp_nondim * eos.R
end

module IAPWS
using parameterized_eos
function gen_iapws()
    h2o_n₀ = [-8.3204464837497, 6.6832105275932, 3.00632, 0.012436, 0.97315, 1.27950,
              0.96956, 0.24873]
    h2o_γ₀ = [1.28728967, 3.53734222, 7.74073708, 9.24437796, 27.5075105]

    h2o_n₁ = [0.12533547935523e-1,
              0.78957634722828e1,
              -0.87803203303561e1,
              0.31802509345418,
              -0.26145533859358,
              -0.78199751687981e-2,
              0.88089493102134e-2]
    h2o_d₁ = [1, 1, 1, 2, 2, 3, 4]
    h2o_t₁ = [-0.5, 0.875, 1, 0.5, 0.75, 0.375, 1]

    h2o_n₂ = [-0.66856572307965,
              0.20433810950965,
              -0.66212605039687e-4,
              -0.19232721156002,
              -0.25709043003438,
              0.16074868486251,
              -0.40092828925807e-1,
              0.39343422603254e-6,
              -0.75941377088144e-5,
              0.56250979351888e-3,
              -0.15608652257135e-4,
              0.11537996422951e-8,
              0.36582165144204e-6,
              -0.13251180074668e-11,
              -0.62639586912454e-9,
              -0.10793600908932,
              0.17611491008752e-1,
              0.22132295167546,
              -0.40247669763528,
              0.58083399985759,
              0.49969146990806e-2,
              -0.31358700712549e-1,
              -0.74315929710341,
              0.47807329915480,
              0.20527940895948e-1,
              -0.13636435110343,
              0.14180634400617e-1,
              0.83326504880713e-2,
              -0.29052336009585e-1,
              0.38615085574206e-1,
              -0.20393486513704e-1,
              -0.16554050063734e-2,
              0.19955571979541e-2,
              0.15870308324157e-3,
              -0.16388568342530e-4,
              0.43613615723811e-1,
              0.34994005463765e-1,
              -0.76788197844621e-1,
              0.22446277332006e-1,
              -0.62689710414685e-4,
              -0.55711118565645e-9,
              -0.19905718354408,
              0.31777497330738,
              -0.11841182425981]
    h2o_c₂ = [1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 2,2,2,2,2,2,2,2, 2,2,2,2,2,2,2,2,2,
              2,2,2, 3,3,3,3, 4,6,6,6,6]
    h2o_d₂ = [1,1,1,2,2,3,4,4,5,7,9,10,11,13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,
              10,10,12,3,4,4,5,14,3,6,6,6]
    h2o_t₂ = [4,6,12,1,5,4,2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,
              4,8,6,9,8,16,22,23,23,10,50,44,46,50]

    h2o_d₃ = [3, 3, 3]
    h2o_t₃ = [0,1,4]
    h2o_n₃ = [-0.31306260323435e2, 0.31546140237781e2, -0.25213154341695e4]
    h2o_α₃ = [20,20,20]
    h2o_β₃ = [150,150,150]
    h2o_γ₃ = [1.21,1.21,1.25]
    h2o_ϵ₃ = [1,1,1]

    h2o_a₄ = [3.5,3.5]
    h2o_b₄ = [0.85,0.95]
    h2o_β₄ = [0.3,0.3]
    h2o_n₄ = [-0.14874640856724, 0.31806110878444]
    h2o_A₄ = [0.32,0.32]
    h2o_B₄ = [0.2,0.2]
    h2o_C₄ = [28,32]
    h2o_D₄ = [700,800]

    Tc = 647.096
    ρc = 322
    R = 0.46151805*1000

    return MultiparameterEOS(Tc,ρc,R,h2o_n₀,h2o_γ₀,h2o_n₁,h2o_d₁,h2o_t₁,
      h2o_n₂,h2o_d₂,h2o_t₂,h2o_c₂,h2o_n₃,h2o_d₃,h2o_t₃,h2o_α₃,h2o_β₃,h2o_γ₃,
      h2o_ϵ₃,h2o_n₄,h2o_a₄,h2o_b₄,h2o_β₄,h2o_A₄,h2o_B₄,h2o_C₄,h2o_D₄)
end
end # ends module IAPWS

module co2
import parameterized_eos
function gen_co2()
    co2_n₀= [ 8.37304456, -3.70454304,  2.5       ,  1.99427042,  0.62105248,
        0.41195293,  1.04028922,  0.08327678]

    co2_γ₀ = [3.15163,   6.1119 ,   6.77708,
            11.32384,  27.08792]

    co2_n₁ = [ 0.38856823,  2.93854759, -5.58671885, -0.767532  ,  0.31729004,
            0.54803316,  0.12279411]
    co2_d₁ = [ 1.,  1.,  1.,  1.,  2.,  2.,  3.]
    co2_t₁ = [ 0.  ,  0.75,  1.  ,  2.  ,  0.75,  2.  ,  0.75]
    co2_n₂ = [  2.16589615e+00,   1.58417351e+00,  -2.31327054e-01,
             5.81169164e-02,  -5.53691372e-01,   4.89466159e-01,
            -2.42757398e-02,   6.24947905e-02,  -1.21758602e-01,
            -3.70556853e-01,  -1.67758797e-02,  -1.19607366e-01,
            -4.56193625e-02,   3.56127893e-02,  -7.44277271e-03,
            -1.73957049e-03,  -2.18101213e-02,   2.43321666e-02,
            -3.47401334e-02,   1.43387158e-01,  -1.34919691e-01,
            -2.31512251e-02,   1.23631255e-02,   2.10583220e-03,
            -3.39585190e-04,   5.59936518e-03,  -3.03351181e-04]

    co2_d₂ = [  1.,   2.,   4.,   5.,   5.,   5.,   6.,   6.,   6.,   1.,   1.,
             4.,   4.,   4.,   7.,   8.,   2.,   3.,   3.,   5.,   5.,   6.,
             7.,   8.,  10.,   4.,   8.]

    co2_c₂ = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  2.,  2.,  2.,  2.,
            2.,  2.,  2.,  3.,  3.,  3.,  4.,  4.,  4.,  4.,  4.,  4.,  5.,  6.]

    co2_t₂ = [  1.5,   1.5,   2.5,   0. ,   1.5,   2. ,   0. ,   1. ,   2. ,
             3. ,   6. ,   3. ,   6. ,   8. ,   6. ,   0. ,   7. ,  12. ,
            16. ,  22. ,  24. ,  16. ,  24. ,   8. ,   2. ,  28. ,  14. ]

    co2_n₃ = [  -213.65488688,  26641.56914927, -24027.21220456,  -283.41603424,
              212.472844  ]
    co2_n₄ = [-0.66642277,  0.72608632,  0.05506867]
    co2_d₃ = [2.,2.,2.,3.,3.]
    co2_t₃ = [1.,0.,1.,3.,3.]
    co2_α₃ = [25.,25.,25.,15.,20.]
    co2_β₃ = [325.,300.,300.,275.,275.]
    co2_γ₃ = [1.16,1.19,1.19,1.25,1.22]
    co2_ϵ₃ = [1.,1.,1.,1.,1.]

    co2_a₄ = [3.5,3.5,3.0]
    co2_b₄ = [0.875,0.925,0.875]
    co2_β₄ = [0.3,0.3,0.3]
    co2_A₄ = [0.7,0.7,0.7]
    co2_B₄ = [0.3,0.3,1.0]
    co2_C₄ = [10.0,10.0,12.5]
    co2_D₄ = [275,275,275]

    Tc = 304.1282
    ρc = 467.6
    pc = 7.38e6
    R = 188.9241
    ω = 0.228
    return MultiparameterEOS(Tc,ρc,pc,ω,R,co2_n₀,co2_γ₀,co2_n₁,co2_d₁,co2_t₁,
      co2_n₂,co2_d₂,co2_t₂,co2_c₂,co2_n₃,co2_d₃,co2_t₃,co2_α₃,co2_β₃,co2_γ₃,
      co2_ϵ₃,co2_n₄,co2_a₄,co2_b₄,co2_β₄,co2_A₄,co2_B₄,co2_C₄,co2_D₄)
end
#eos = gen_co2()
#Tc = 304.1282
#ρc = 467.6
#δ₁ = 838.025/ρc
#τ₁ = Tc/500
#test_peos(eos, δ₁, τ₁)
end #ends module co2


text:
using RandomMatrices
using Base.Test

let
    d = Semicircle()
    @test cdf(d, 2.0) == 1.0
    @test cdf(d, -2.0) == 0.0
    @test pdf(d, 2.0) == 0.0
    @test pdf(d, -2.0) == 0.0

    @test entropy(d) == log(2π) - 0.5
    @test mean(d) == median(d) == skewness(d) == 0
    @test modes(d) == [0]
    @test std(d) == var(d) == 1
    @test kurtosis(d) == 2

    @test isfinite(rand(d))
    @test moment(d, 4) == 2
    @test cumulant(d, 4) == 1
    @test moment(d, 5) == cumulant(d, 5) == freecumulant(d, 5) == 0

    @test cf(d, 1.0) == 0.0

    #XXX Throws AmosException - ???
    #@show mgf(d, 1.0)

    x = rand(d)
    @test insupport(d, x)
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
using HDF5

if length(ARGS) >= 1
    batch_size = parse(Int, ARGS[1])
else
    batch_size = 64
end
net = Net(batch_size)
data, label  = HDF5DataLayer(net, "data/train.txt", "data/val.txt")
trans = TransformLayer(net, :trans, data; crop=(224, 224), random_mirror=true)

conv1_1 = ConvolutionLayer( :conv1_1, net, trans,   64,  3, 1, 1)
relu1_1 = ReLULayer(        :relu1_1, net, conv1_1)
pool1   = MaxPoolingLayer(  :pool1,   net, relu1_1, 2,   2, 0)

conv2_1 = ConvolutionLayer( :conv2_1, net, pool1,   128, 3, 1, 1)
relu2_1 = ReLULayer(        :relu2_1, net, conv2_1)
pool2   = MaxPoolingLayer(  :pool2,   net, relu2_1, 2,   2, 0)

conv3_1 = ConvolutionLayer( :conv3_1, net, pool2,   256, 3, 1, 1)
relu3_1 = ReLULayer(        :relu3_1, net, conv3_1)
conv3_2 = ConvolutionLayer( :conv3_2, net, relu3_1, 256, 3, 1, 1)
relu3_2 = ReLULayer(        :relu3_2, net, conv3_2)
pool3   = MaxPoolingLayer(  :pool3,   net, relu3_2, 2,   2, 0)

conv4_1 = ConvolutionLayer( :conv4_1, net, pool3,   512, 3, 1, 1)
relu4_1 = ReLULayer(        :relu4_1, net, conv4_1)
conv4_2 = ConvolutionLayer( :conv4_2, net, relu4_1, 512, 3, 1, 1)
relu4_2 = ReLULayer(        :relu4_2, net, conv4_2)
pool4   = MaxPoolingLayer(  :pool4,   net, relu4_2, 2,   2, 0)

conv5_1 = ConvolutionLayer( :conv5_1, net, pool4,   512, 3, 1, 1)
relu5_1 = ReLULayer(        :relu5_1, net, conv5_1)
conv5_2 = ConvolutionLayer( :conv5_2, net, relu5_1, 512, 3, 1, 1)
relu5_2 = ReLULayer(        :relu5_2, net, conv5_2)
pool5   = MaxPoolingLayer(  :pool5,   net, relu5_2, 2,   2, 0)

fc6     = InnerProductLayer(:fc6,     net, pool5, 4096)
relu6   = ReLULayer(        :relu6,   net, fc6)
drop6   = DropoutLayer(     :drop6,   net, relu6, .5f0)
fc7     = InnerProductLayer(:fc7,     net, drop6, 4096)
relu7   = ReLULayer(        :relu7,   net, fc7)
drop7   = DropoutLayer(     :drop7,   net, relu7, .5f0)
fc8     = InnerProductLayer(:fc8,     net, drop7, 1000)

loss     = SoftmaxLossLayer(:loss, net, fc8, label)
accuracy = AccuracyLayer(:accuracy, net, fc8, label)

params = SolverParameters(
    lr_policy    = LRPolicy.Decay(.01f0, 5.0f-7),
    mom_policy   = MomPolicy.Fixed(0.9),
    max_epoch    = 300,
    regu_coef    = .0005)
sgd = SGD(params)
solve(sgd, net)


text:
using MCMCTempering
using Test
using Distributions
using Plots
using AdvancedMH
using MCMCChains

@testset "MCMCTempering.jl" begin

    # θᵣ = [-1., 1., 2., 1., 15., 2., 90., 1.5]
    # γs = [0.15, 0.25, 0.3, 0.3]

    # Δ = check_Δ([0, 0.01, 0.1, 0.25, 0.5, 1])

    # modelᵣ = MixtureModel(Distributions.Normal.(eachrow(reshape(θᵣ, (2,4)))...), γs)
    # # xrange = -10:0.1:100
    # # tempered_densities = pdf.(modelᵣ, xrange) .^ Δ'
    # # norm_const = sum(tempered_densities[:,1])
    # # for i in 2:length(Δ)
    # #     tempered_densities[:,i] = (tempered_densities[:,i] ./ sum(tempered_densities[:,i])) .* norm_const
    # # end
    # # plot(xrange, tempered_densities, label = Δ')

    # data = rand(modelᵣ, 100)

    # insupport(θ) = all(reshape(θ, (2,4))[2,:] .≥ 0)
    # dist(θ) = MixtureModel(Distributions.Normal.(eachrow(reshape(θ, (2,4)))...), γs)
    # density(θ) = insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf

    # # Construct a DensityModel.
    # model = DensityModel(density)

    # # Set up our sampler with a joint multivariate Normal proposal.
    # spl = RWMH(MvNormal(8,1))

    # @test chain, temps = SimulatedTempering(model, spl, Δ, chain_type=Chains)
    # @test chains, temps = ParallelTempering(model, spl, Δ, chain_type=Chains)

end