text:
struct ConstraintKey{T <: ConstraintType, U <: Union{PSY.Component, PSY.System}} <:
       OptimizationContainerKey
    meta::String
end

function ConstraintKey(
    ::Type{T},
    ::Type{U},
    meta=CONTAINER_KEY_EMPTY_META,
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}}
    check_meta_chars(meta)
    return ConstraintKey{T, U}(meta)
end

get_entry_type(
    ::ConstraintKey{T, U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}} = T
get_component_type(
    ::ConstraintKey{T, U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}} = U

function encode_key(key::ConstraintKey)
    return encode_symbol(get_component_type(key), get_entry_type(key), key.meta)
end

Base.convert(::Type{ConstraintKey}, name::Symbol) = ConstraintKey(decode_symbol(name)...)

struct AbsoluteValueConstraint <: ConstraintType end
struct ActiveConstraint <: ConstraintType end
struct ActiveRangeConstraint <: ConstraintType end
struct ActiveRangeICConstraint <: ConstraintType end
struct AreaDispatchBalanceConstraint <: ConstraintType end
struct AreaParticipationAssignmentConstraint <: ConstraintType end
struct BalanceAuxConstraint <: ConstraintType end
struct CommitmentConstraint <: ConstraintType end
struct CopperPlateBalanceConstraint <: ConstraintType end
struct DurationConstraint <: ConstraintType end
struct EnergyBalanceConstraint <: ConstraintType end
struct EnergyBudgetConstraint <: ConstraintType end
struct EnergyCapacityConstraint <: ConstraintType end
struct EnergyCapacityDownConstraint <: ConstraintType end
struct EnergyCapacityUpConstraint <: ConstraintType end
struct EnergyLimitConstraint <: ConstraintType end
struct EnergyTargetConstraint <: ConstraintType end
struct EnergyShortageVariableLimitsConstraint <: ConstraintType end
struct EqualityConstraint <: ConstraintType end
struct FeedforwardSemiContinousConstraint <: ConstraintType end
struct FeedforwardIntegralLimitConstraint <: ConstraintType end
struct FeedforwardUpperBoundConstraint <: ConstraintType end
struct FeedforwardLowerBoundConstraint <: ConstraintType end
struct FeedforwardEnergyTargetConstraint <: ConstraintType end
struct FlowActivePowerConstraint <: ConstraintType end
struct FlowActivePowerFromToConstraint <: ConstraintType end
struct FlowActivePowerToFromConstraint <: ConstraintType end
struct FlowLimitConstraint <: ConstraintType end
struct FlowLimitFromToConstraint <: ConstraintType end
struct FlowLimitToFromConstraint <: ConstraintType end
struct FlowRateConstraint <: ConstraintType end
struct FlowRateConstraintFromTo <: ConstraintType end
struct FlowRateConstraintToFrom <: ConstraintType end
struct FlowReactivePowerConstraint <: ConstraintType end
struct FlowReactivePowerFromToConstraint <: ConstraintType end
struct FlowReactivePowerToFromConstraint <: ConstraintType end
struct HVDCPowerBalance <: ConstraintType end
struct FrequencyResponseConstraint <: ConstraintType end
struct InflowRangeConstraint <: ConstraintType end
struct InputPowerRangeConstraint <: ConstraintType end
struct InterConnectionLimitConstraint <: ConstraintType end
struct MustRunConstraint <: ConstraintType end
struct NetworkFlowConstraint <: ConstraintType end
struct NodalBalanceActiveConstraint <: ConstraintType end
struct NodalBalanceReactiveConstraint <: ConstraintType end
struct PowerOutputRangeConstraint <: ConstraintType end
struct ParticipationAssignmentConstraint <: ConstraintType end
struct PieceWiseLinearCostConstraint <: ConstraintType end
struct RampConstraint <: ConstraintType end
struct RampLimitConstraint <: ConstraintType end
struct RangeLimitConstraint <: ConstraintType end
struct RateLimitConstraint <: ConstraintType end
struct RateLimitConstraintFromTo <: ConstraintType end
struct RateLimitConstraintToFrom <: ConstraintType end
struct ReactivePowerConstraint <: ConstraintType end
struct ReactiveRangeConstraint <: ConstraintType end
struct RegulationLimitsConstraint <: ConstraintType end
struct RequirementConstraint <: ConstraintType end
struct ReserveEnergyConstraint <: ConstraintType end
struct ReservePowerConstraint <: ConstraintType end
struct ComponentReserveUpBalance <: ConstraintType end
struct ComponentReserveDownBalance <: ConstraintType end
struct SACEPIDAreaConstraint <: ConstraintType end
struct StartTypeConstraint <: ConstraintType end
struct StartupInitialConditionConstraint <: ConstraintType end
struct StartupTimeLimitTemperatureConstraint <: ConstraintType end

abstract type PowerVariableLimitsConstraint <: ConstraintType end
struct EnergyVariableLimitUpConstraint <: ConstraintType end
struct EnergyVariableLimitDownConstraint <: ConstraintType end
struct InputActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end
struct OutputActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end
struct ActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end
struct ReactivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end
struct ActivePowerVariableTimeSeriesLimitsConstraint <: PowerVariableLimitsConstraint end
struct ComponentActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end
struct ComponentReactivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end

# These apply to the processing of constraint duals
should_write_resulting_value(::Type{<:ConstraintType}) = true
convert_result_to_natural_units(::Type{<:ConstraintType}) = false


text:
#using PyPlot
#addprocs(7)
using ProgressMeter
@everywhere using DataFrames

@everywhere type State
    current_state::Int64
    policy::Vector{Float64}
    nr_world_states::Int64
    nr_actions::Int64
    target_sequence::Vector{Int64}
    sum_of_rewards::Vector{Float64}
    nr_of_evaluations::Vector{Float64}
    value_function::Vector{Float64} # on the world state
    state_action_value_function::Matrix{Float64} # on the world state
    k::Int64
    value_history::Matrix{Float64}
end

@everywhere type Result
    experiment::State
    values::Matrix{Float64}
end

@everywhere function create_experiment(nr_world_states::Int64, nr_actions::Int64, target_actions::Int64, k::Int64)
    policy = rand(nr_actions)
    policy = policy / sum(policy)

    min = minimum([target_actions, nr_world_states])
    max = maximum([0, nr_world_states - nr_actions])

    sequence = vcat([1:min], int(ceil(rand(max) * nr_actions)))
    sequence = sequence[randperm(length(sequence))]

    return State(1, policy, nr_world_states, nr_actions, sequence, zeros(nr_world_states), zeros(nr_world_states), zeros(nr_world_states), zeros(nr_world_states, nr_actions), k, zeros(2,2))
end

@everywhere function run_episode(state::State, initial_state::Int64, T::Int64)
    state.current_state = initial_state
    for t = 1:T
        p = rand()
        s = 0
        action = state.nr_actions
        for i = 1:state.nr_actions
            s = s + state.policy[i]
            if p < s
                action = i
                break
            end
        end
        if action == state.target_sequence[state.current_state]
            state.current_state = state.current_state + 1
        else
            state.current_state = 1
        end
        if state.current_state >= state.nr_world_states
            return 1
        end
    end
    return 0
end

@everywhere function monte_carlo_estimation_of_world_state_function!(state::State, T::Int64, nr_of_episodes::Int64)
    for e = 1:nr_of_episodes
        initial_state                          = int64(ceil(rand()*state.nr_world_states))
        if initial_state == 1
          println("initial state: $(initial_state)")
        end
        reward                                 = run_episode(state, initial_state, T)
        state.sum_of_rewards[initial_state]    = state.sum_of_rewards[initial_state] + reward
        state.nr_of_evaluations[initial_state] = state.nr_of_evaluations[initial_state] + 1.0
        state.value_function[initial_state]    = state.sum_of_rewards[initial_state] / state.nr_of_evaluations[initial_state]
    end
end

@everywhere function calculate_world_state_action_function!(state::State)
    for i = 1:state.nr_world_states
        for action = 1:state.nr_actions
            if action == state.target_sequence[i]
                if i == state.nr_world_states
                    state.state_action_value_function[i,action] = 1.0
                else
                    state.state_action_value_function[i,action] = state.value_function[i+1]
                end
            else
                state.state_action_value_function[i,action] = state.value_function[1]
            end
        end
    end

    for action = 1:state.nr_actions
        if action == state.target_sequence[state.nr_world_states]
            state.state_action_value_function[state.nr_world_states,action] = 1.0
        else
            state.state_action_value_function[state.nr_world_states,action] = state.value_function[1]
        end
    end
end

@everywhere function prune_world_state_action_function!(state::State)
    nr_of_values = minimum([state.k, state.nr_actions])
    for i = 1:state.nr_world_states
        if sum(state.state_action_value_function[i,:]) > 0.0
            values = DataFrame(VALUES=[v for v in state.state_action_value_function[i,:]], INDICES=[1:state.nr_actions])
            sort!(values, cols = (:VALUES), rev=true)

            values = values[1:nr_of_values,:]
            for a = 1:state.nr_actions
                state.state_action_value_function[i,a] = 0.0
            end

            for a = 1:size(values)[1]
                state.state_action_value_function[i,values[:INDICES][a]] = values[:VALUES][a]
            end
        else
            indices = randperm(state.nr_actions)[1:nr_of_values]
            for a in indices
                state.state_action_value_function[i, a] = 0.001
            end
        end
    end
end

@everywhere function update_policy_from_world_state_action_function!(state::State)
    sum_of_values = zeros(state.nr_actions)
    for a = 1:state.nr_actions
        sum_of_values[a] = sum(state.state_action_value_function[:,a])
    end

    s = sum(sum_of_values)
    for a = 1:state.nr_actions
        state.policy[a] = sum_of_values[a] / s
    end
end

@everywhere function update_policy!(state::State)
    calculate_world_state_action_function!(state)
    prune_world_state_action_function!(state)
    update_policy_from_world_state_action_function!(state)
end

@everywhere function scan_over_k(k::Int64, N::Int64, episode_length::Int64, nr_of_episodes::Int64)
    nr_world_states   = k + 1
    nr_of_actions     = k
    exp               = create_experiment(nr_world_states, nr_of_actions, nr_of_actions, k)
    exp.value_history = zeros(N, nr_world_states)
    pm = Progress(N, 1)
    for i = 1:N
        monte_carlo_estimation_of_world_state_function!(exp, (episode_length+1) * k, nr_of_episodes)
        update_policy!(exp)
        exp.value_history[i,:] = exp.value_function
        next!(pm)
    end
    return exp
end

#= N               = 1000000 =#
#= nr_world_states = 5 =#
#= nr_of_actions   = 5 =#
#= k               = 5 =#
#= v,r = scan_over_k(k, N, nr_world_states, nr_of_actions) =#

#= println("policy") =#
#= println(r.policy) =#

#= println("state action value function") =#
#= println(r.state_action_value_function) =#

ext: jl

max_stars_repo_name: kzahedi/POMDP

max_issues_repo_name: kzahedi/POMDP

max_forks_repo_name: kzahedi/POMDP


text:
using JToric
using Test

@testset "JToric.jl" begin
    # load necessary gap-packages
    using CapAndHomalg
    GAP.Globals.LoadPackage( GapObj( "JConvex" ) )
    CapAndHomalg.LoadPackage( "ToricVarieties" )

    # Compute properties of toric varieties on the example of a Hirzebruch surface
    Rays = [[-1,5],[0,1],[1,0],[0,-1]]
    Cones = [[1,2],[2,3],[3,4],[4,1]]
    H5 = toric_variety( Rays, Cones )
    @test JToric.is_normal_variety( H5 ) == true
    @test JToric.is_affine( H5 ) == false
    @test JToric.is_projective( H5 ) == true
    @test JToric.is_smooth( H5 ) == true
    @test JToric.is_complete( H5 ) == true
    @test JToric.has_torusfactor( H5 ) == false
    @test JToric.has_no_torusfactor( H5 ) == true
    @test JToric.is_orbifold( H5 ) == true
    @test JToric.is_simplicial( H5 ) == true
    @test JToric.is_isomorphic_to_projective_space( H5 ) == false
    @test JToric.is_direct_product_of_projective_spaces( H5 ) == false

    # Compute properties of toric varieties on the example of a projective space
    P2 = JToric.projective_space( 2 )
    @test JToric.is_normal_variety( P2 ) == true
    @test JToric.is_affine( P2 ) == false
    @test JToric.is_projective( P2 ) == true
    @test JToric.is_smooth( P2 ) == true
    @test JToric.is_complete( P2 ) == true
    @test JToric.has_torusfactor( P2 ) == false
    @test JToric.has_no_torusfactor( P2 ) == true
    @test JToric.is_orbifold( P2 ) == true
    @test JToric.is_simplicial( P2 ) == true
    @test JToric.is_isomorphic_to_projective_space( P2 ) == true
    @test JToric.is_direct_product_of_projective_spaces( P2 ) == true

    # compute properties of toric divisors on the example of the trivial divisor
    D=create_divisor( [ 0,0,0,0 ], H5 )
    @test JToric.is_cartier( D ) == true
    @test JToric.is_principal( D ) == true
    @test JToric.is_primedivisor( D ) == false
    @test JToric.is_basepoint_free( D ) == true
    @test JToric.is_ample( D ) == false
    @test JToric.is_very_ample( D ) == "fail"
    @test JToric.is_numerically_effective( D ) == true

    # compute properties of toric divisors on the example of a non-trivial divisor
    D2 = divisor_of_character( [ 1,2 ], H5 )
    @test JToric.is_cartier( D2 ) == true
    @test JToric.is_principal( D2 ) == true
    @test JToric.is_primedivisor( D2 ) == false
    @test JToric.is_basepoint_free( D2 ) == true
    @test JToric.is_ample( D2 ) == false
    @test JToric.is_very_ample( D2 ) == "fail"
    @test JToric.is_numerically_effective( D2 ) == true

    # compute properties of toric divisors on the example of another non-trivial divisor
    D3 = divisor_of_class( H5, [ 1,2 ] )
    @test JToric.is_cartier( D3 ) == true
    @test JToric.is_principal( D3 ) == false
    @test JToric.is_primedivisor( D3 ) == false
    @test JToric.is_basepoint_free( D3 ) == true
    @test JToric.is_ample( D3 ) == true
    @test JToric.is_very_ample( D3 ) == true
    @test JToric.is_numerically_effective( D3 ) == true

    # compute attributes of toric varieties on the example of the Hirzebruch surface H5 defined above
    cover = affine_open_covering( H5 )
    @test size( cover )[ 1 ] == 4
    cox_ring( H5 )
    @test size( list_of_variables_of_cox_ring( H5 ) )[ 1 ] == 4
    class_group( H5 )
    torus_invariant_divisor_group( H5 )
    map_from_character_to_principal_divisor( H5 )
    map_from_weil_divisors_to_class_group( H5 )
    @test dimension( H5 ) == 2
    @test dimension_of_torusfactor( H5 ) == 0
    coordinate_ring_of_torus( H5 )
    @test size( list_of_variables_of_coordinate_ring_of_torus( H5 ) )[ 1 ] == 4
    @test size( is_product_of( H5 ) )[ 1 ] == 1
    character_lattice( H5 )
    divisors = torus_invariant_prime_divisors( H5 )
    @test ( 0 in [ is_primedivisor( d ) for d in divisors ] ) == false
    irrelevant_ideal( H5 )
    stanley_reisner_ideal( H5 )
    morphism_from_cox_variety( H5 )
    cox_variety( H5 )
    fan_of_variety( H5 )
    fan( H5 )
    cartier_torus_invariant_divisor_group( H5 )
    picard_group( H5 )
    @test name_of_variety( H5 ) == "No name set for this variety"
    set_name_of_variety( H5, "Hirzebruch surface" )
    @test name_of_variety( H5 ) == "Hirzebruch surface"
    zariski_cotangent_sheaf( H5 )
    cotangent_sheaf( H5 )
    @test euler_characteristic( H5 ) == 0
    # UnderlyingSheaf( H5 ) <- Error in gap

    # apply methods to toric varieties on the example of the Hirzebruch surface H5 and projective space P2 defined above
    coordinate_ring_of_torus( H5, [ "u", "v", "w", "z" ] )
    cox_ring( H5, "u" )
    v = H5 * P2;
    @test JToric.is_normal_variety( v ) == true
    @test JToric.is_affine( v ) == false
    @test JToric.is_projective( v ) == true
    @test JToric.is_smooth( v ) == true
    @test JToric.is_complete( v ) == true
    @test JToric.has_torusfactor( v ) == false
    @test JToric.has_no_torusfactor( v ) == true
    @test JToric.is_orbifold( v ) == true
    @test JToric.is_simplicial( v ) == true
    @test JToric.is_isomorphic_to_projective_space( v ) == false
    @test JToric.is_direct_product_of_projective_spaces( v ) == false
end

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
using Documenter, TupleTools

makedocs(modules=[TupleTools],
            sitename="TupleTools.jl",
            authors = "Jutho Haegeman",
            format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
            pages = [
                "Home" => "index.md",
            ])

deploydocs(repo = "github.com/Jutho/TupleTools.jl.git")

max_stars_repo_name: Jutho/TupleTools

max_issues_repo_name: Jutho/TupleTools

max_forks_repo_name: Jutho/TupleTools


text:
fsa = true

using StaticArrays
@static if fsa
    using FixedSizeArrays
end

if !isdefined(:N)
    N = 4
end

M_f = div(10^9, N^3)
M_g = div(2*10^8, N^2)

# Size

A = rand(N,N)
As = SMatrix{N,N}(A)
Am = MMatrix{N,N}(A)
@static if fsa
    Af = Mat(ntuple(j -> ntuple(i->A[i,j], N), N)) # there is a bug in FixedSizeArrays Mat constructor (13 July 2016)
end

if !isdefined(:f_mut_marray) || !isdefined(:benchmark_suite) || benchmark_suite == false
    f(n::Integer, A) = @inbounds (C = A; for i = 1:n; C = C*A; end; return C)
    f_unrolled{M}(n::Integer, A::Union{SMatrix{M,M},MMatrix{M,M}}) = @inbounds (C = A; for i = 1:n; C = StaticArrays.A_mul_B_unrolled(C,A); end; return C)
    f_unrolled_chunks{M}(n::Integer, A::Union{SMatrix{M,M},MMatrix{M,M}}) = @inbounds (C = A; for i = 1:n; C = StaticArrays.A_mul_B_unrolled_chunks(C,A); end; return C)
    f_via_sarray{M}(n::Integer, A::MMatrix{M,M})= @inbounds (C = A; for i = 1:n; C = MMatrix{M,M}(SMatrix{M,M}(C)*SMatrix{M,M}(A)); end; return C)
    f_mut_array(n::Integer, A) = @inbounds (C = copy(A); tmp = similar(A); for i = 1:n;  A_mul_B!(tmp, C, A); map!(identity, C, tmp); end; return C)
    f_mut_marray(n::Integer, A) = @inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B_unrolled!(tmp, C, A); C.data = tmp.data; end; return C)
    f_blas_marray(n::Integer, A) = @inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B_blas!(tmp, C, A); C.data = tmp.data; end; return C)

    g(n::Integer, A) = @inbounds (C = A; for i = 1:n; C = C + A; end; return C)
    g_mut(n::Integer, A) = @inbounds (C = similar(A); C[:] = A[:]; for i = 1:n; @inbounds map!(+, C, C, A); end; return C)
    g_via_sarray{M}(n::Integer, A::MMatrix{M,M}) = @inbounds (C = similar(A); C[:] = A[:]; for i = 1:n; C = MMatrix{M,M}(SMatrix{M,M}(C) + SMatrix{M,M}(A)); end; return C)

    # Notes: - A[:] = B[:] is allocating in Base, unlike `map!`
    #        - Also, the same goes for Base's implementation of broadcast!(f, A, B, C) (but map! is OK).
    #        - I really need to implement copy() in StaticArrays... (and maybe a special instance of setindex!(C, :))

    # Old, "broken" defs:
    #f(n::Integer, A, B) = @inbounds (C = A*B; for i = 1:n; C += A*B; end; return C)
    #f_mut(n::Integer, A, B) = @inbounds (C = A*B; tmp = similar(C); for i = 1:n; A_mul_B!(tmp, A, B); map!(+, C, C, tmp); end; return C)
    #g(n::Integer, A, B) = @inbounds (C = A + B; for i = 1:n; C = C + B; end; return C)
    #g_mut(n::Integer, A, B) = @inbounds (C = A + B; for i = 1:n; @inbounds broadcast!(+, C, C, B); end; return C)
end

# Don't count these compilations either
C = f(2, A)
C_mut = f_mut_array(2, A)
if N == 4;  @assert C_mut ≈ C; end

C = g(2, A)
C_mut = g_mut(2, A)

if N == 4; @assert C_mut ≈ C; end

println("=====================================")
println("    Benchmarks for $N×$N matrices")
println("=====================================")
print("StaticArrays compilation time (×3):")
@time eval(quote
    # Warmup and some checks
    C = f(2, A)
    Cs = f(2, As)
    Cs_unrolled = f_unrolled(2, As)
    Cs_chunks = f_unrolled_chunks(2, As)
    Cm = f(2, Am)
    Cm_unrolled = f_unrolled(2, Am)
    Cm_chunks = f_unrolled_chunks(2, Am)
    Cm_via_sarray = f_via_sarray(2, Am)
    Cm_mut = f_mut_marray(2, Am)
    Cm_blas = f_blas_marray(2, Am)

    Cs::SMatrix
    if N == 4; @assert Cs ≈ C; end
    Cs_unrolled::SMatrix
    if N == 4; @assert Cs_unrolled ≈ C; end
    Cs_chunks::SMatrix
    if N == 4; @assert Cs_chunks ≈ C; end
    Cm::MMatrix
    if N == 4; @assert Cm ≈ C; end
    Cm_unrolled::MMatrix
    if N == 4; @assert Cm_unrolled ≈ C; end
    Cm_chunks::MMatrix
    if N == 4; @assert Cm_chunks ≈ C; end
    Cm_via_sarray::MMatrix
    if N == 4; @assert Cm_via_sarray ≈ C; end
    Cm_mut::MMatrix
    if N == 4; @assert Cm_mut ≈ C; end
    Cm_blas::MMatrix
    if N == 4; @assert Cm_blas ≈ C; end

    C = g(2, A)
    Cs = g(2, As)
    Cm = g(2, Am)
    Cm_via_sarray = g_via_sarray(2, Am)
    Cm_mut = g_mut(2, Am)

    Cs::SMatrix
    if N == 4; @assert Cs == C; end
    Cm::MMatrix
    if N == 4; @assert Cm == C; end
    Cm_via_sarray::MMatrix
    if N == 4; @assert Cm_via_sarray == C; end
    Cm_mut::MMatrix
    if N == 4; @assert Cm_mut == C; end
end)
@static if fsa
    print("FixedSizeArrays compilation time:  ")
    @time eval(quote
        C = f(2, A)
        Cf = f(2, Af)
        Cf::Mat
        if N == 4; @assert Cf == C; end

        C = g(2, A)
        Cf = g(2, Af)
        Cf::Mat
        if N == 4; @assert Cf == C; end
    end)
end
println()

# Do the performance tests

println("Matrix multiplication")
println("---------------------")
begin
   print("Array               ->"); @time f(M_f, A)
   print("Array (mutating)    ->"); @time f_mut_array(M_f, A)
   print("SArray              ->"); @time f(M_f, As)
   print("SArray (unrolled)   ->"); @time f_unrolled(M_f, As)
   print("SArray (chunks)     ->"); @time f_unrolled_chunks(M_f, As)
   print("MArray              ->"); @time f(M_f, Am)
   print("MArray (unrolled)   ->"); @time f_unrolled(M_f, Am)
   print("MArray (chunks)     ->"); @time f_unrolled_chunks(M_f, Am)
   print("MArray (via SArray) ->"); @time f_via_sarray(M_f, Am)
   print("MArray (mutating)   ->"); @time f_mut_marray(M_f, Am)
   print("MArray (BLAS gemm!) ->"); @time f_blas_marray(M_f, Am)
   @static if fsa
       print("Mat                 ->"); @time f(M_f, Af)
   end
end
println()

println("Matrix addition")
println("---------------")
begin
   print("Array               ->"); @time g(M_g, A)
   print("Array (mutating)    ->"); @time g_mut(M_g, A) # broadcast! seems to be broken!
   print("SArray (unrolled)   ->"); @time g(M_g, As)
   print("MArray (unrolled)   ->"); @time g(M_g, Am)
   print("MArray (via SArray) ->"); @time g_via_sarray(M_g, Am)
   print("MArray (mutating)   ->"); @time g_mut(M_g, Am)
   @static if fsa
       print("Mat                 ->"); @time g(M_g, Af)
   end
end
println()


text:
function gravityFridaLS2(q)
    baseAngles  = [-0.63 , 0.95 , -0.18]
    g = 9.82
    c1 = cos(q[1]);
    s1 = sin(q[1]);
    c2 = cos(q[2]);
    s2 = sin(q[2]);
    c3 = cos(q[3]);
    s3 = sin(q[3]);
    c4 = cos(q[4]);
    s4 = sin(q[4]);
    c5 = cos(q[5]);
    s5 = sin(q[5]);
    c6 = cos(q[6]);
    s6 = sin(q[6]);
    c7 = cos(q[7]);
    s7 = sin(q[7]);

    sb1 = sin(baseAngles[1]);
    cb1 = cos(baseAngles[1]);
    sb2 = sin(baseAngles[2]);
    cb2 = cos(baseAngles[2]);
    sb3 = sin(baseAngles[3]);
    cb3 = cos(baseAngles[3]);

    shoulderLen =   0.110;
    shoulderOffs = -0.030;
    elbowOffs =     0.0405;
    wristOffs =     0.0135;
    handOffs =     -0.027;
    upArmLen =      0.2465;
    lowArmLen =     0.265;
    handLen =       0.032;

    A = zeros(7,28);

    p0=g*cb1;
    p1=g*sb1;
    p2=sb2*p0;
    p3=s1*cb3;
    p4=c2*c3;
    p5=cb2*p0;
    p6=c1*cb3;
    p7=c1*sb3;
    p8=s1*sb3;
    p9=s4*p1;
    p10=s4*p2;
    p11=c5*c6;
    p12=c4*p1;
    p13=s1*sb3;
    p14=c4*p2;
    p15=s3*c2;
    p16=c7^2;
    p17=s7^2;
    p18=s5*c6;
    p19=s2*s6;
    p20=s3*s6;
    p21=s2*c3;
    p22=c6^2;
    p23=s6*p4;
    p24=p3*p9;
    p25=p6*p10;
    p26=s4*p5;
    p27=c3*p1;
    p28=p7*p9;
    p29=p8*p10;
    p30=p3*p12;
    p31=p6*p14;
    p32=c4*p5;
    p33=c3*p2;
    p34=p1*p15;
    p35=s6^2;
    p36=p7*p12;
    p37=p8*p14;
    p38=c5*wristOffs;
    p39=p2*p15;
    p40=s2*s3;
    p41=p3*p10;
    p42=p6*p9;
    p43=p7*p10;
    p44=c7*p11;
    p45=s3*s5;
    p46=p3*p14;
    p47=c5*p19;
    p48=p5*p40;
    p49=p9*p13;
    p50=s7*p11;
    p51=p6*p12;
    p52=p7*p14;
    p53=s5*s6;
    p54=p1*p3;
    p55=c3*p19;
    p56=p8*p12;
    p57=handOffs*p11;
    p58=s5*p4;
    p59=s6*c2;
    p60=p1*p7;
    p61=p2*p6;
    p62=p3*p33;
    p63=p2*p8;
    p64=s2*s5;
    p65=c5*c7;
    p66=elbowOffs*p4;
    p67=p6*p27;
    p68=p7*p33;
    p69=c2*p32;
    p70=s7*c5;
    p71=c7*p18;
    p72=s2*c6;
    p73=s3*elbowOffs;
    p74=p8*p27;
    p75=p1*p6;
    p76=p2*p3;
    p77=s5*wristOffs;
    p78=c2*shoulderOffs;
    p79=s7*p18;
    p80=c6*p4;
    p81=p11*p17;
    p82=p1*p8;
    p83=p2*p7;
    p84=p3*p34;
    p85=p6*p39;
    p86=c5*p23;
    p87=elbowOffs*p21;
    p88=s5*c7;
    p89=wristOffs*p45;
    p90=handOffs*p18;
    p91=s6*c6;
    p92=p11*p16;
    p93=p7*p34;
    p94=p8*p39;
    p95=c5*p20;
    p96=s5*s7;
    p97=p21*p26;
    p98=s2*elbowOffs;
    p99=wristOffs*p58;
    p100=p3*p27;
    p101=s3*c6;
    p102=s6*c5;
    p103=p16*p18;
    p104=c2*p26;
    p105=handOffs*p35;
    p106=s2*p30;
    p107=p21*p32;
    p108=p7*p27;
    p109=p6*p34;
    p110=p3*p39;
    p111=p17*p18;
    p112=c5*p4;
    p113=s2*p36;
    p114=p8*p34;
    p115=p6*p33;
    p116=p7*p39;
    p117=p26*p40;
    p118=p32*p55;
    p119=s3*c5;
    p120=s2*p31;
    p121=p4*p24;
    p122=p16*p20;
    p123=p8*p33;
    p124=p40*p54;
    p125=s2*p37;
    p126=p4*p25;
    p127=s3*p42;
    p128=p17*p20;
    p129=p16*p19;
    p130=p15*p24;
    p131=c5*p59;
    p132=p40*p60;
    p133=s2*upArmLen;
    p134=s3*p49;
    p135=p4*p28;
    p136=p17*p19;
    p137=c3*p42;
    p138=s5*p20;
    p139=p40*p61;
    p140=p26*p59;
    p141=s2*shoulderOffs;
    p142=s3*p41;
    p143=p4*p29;
    p144=p17*p23;
    p145=c3*p49;
    p146=p15*p28;
    p147=handOffs*p19;
    p148=c2*upArmLen;
    p149=c6*p11;
    p150=p40*p63;
    p151=p26*p55;
    p152=s3*p43;
    p153=p16*p23;
    p154=s2*p38;
    p155=c3*p41;
    p156=p15*p25;
    p157=c2*p30;
    p158=handOffs*p22;
    p159=c6*p57;
    p160=p17*p50;
    p161=c7*p64;
    p162=p4*p38;
    p163=c3*p43;
    p164=p15*p29;
    p165=s5*p23;
    p166=c2*p36;
    p167=p22*p92;
    p168=p16*p44;
    p169=s7*p64;
    p170=c7*p45;
    p171=s2*lowArmLen;
    p172=s3*p38;
    p173=c3*p47;
    p174=c2*p31;
    p175=p22*p81;
    p176=p16*p50;
    p177=p17*p44;
    p178=handOffs*p23;
    p179=s7*p45;
    p180=c7*p58;
    p181=lowArmLen*p4;
    p182=s6*c3;
    p183=c2*p37;
    p184=s7*p58;
    p185=handOffs*p20;
    p186=p19*p45;
    p187=s6*p15;
    p188=s7*p19;
    p189=p5*p15;
    p190=s6*p92;
    p191=c7*p19;
    p192=p5*p21;
    p193=s6*p81;
    p194=p21*p24;
    p195=p21*p28;
    p196=p21*p25;
    p197=p21*p29;
    A[1,1] = -shoulderOffs*p75+shoulderOffs*p82-shoulderOffs*p76-shoulderOffs*p83;
    A[1,2] = -shoulderOffs*p75+shoulderOffs*p82+p78*p75-p78*p82-shoulderOffs*p76-shoulderOffs*p83+p78*p76+p78*p83;
    A[1,3] = -shoulderOffs*p75+shoulderOffs*p82+p78*p75+p54*p73+p60*p73-p75*p133-p78*p82+p82*p133-shoulderOffs*p76-shoulderOffs*p83+p66*p82-p61*p73+p78*p76+p78*p83+p63*p73-p76*p133-p83*p133-p66*p75-p66*p76-p66*p83;
    A[1,4] = -shoulderOffs*p75+shoulderOffs*p82+p78*p75+p54*p73+p60*p73-p75*p133-p78*p82+p82*p133-shoulderOffs*p76-shoulderOffs*p83+p66*p82+p24*p73+p56*p98+p28*p73-p61*p73+p78*p76+p78*p83+p63*p73-p76*p133-p83*p133-p66*p75-p51*p98-p66*p76-p66*p83-p46*p98-p25*p73-p52*p98+p29*p73-p42*p66+p49*p66-p41*p66-p43*p66;
    A[1,5] = -shoulderOffs*p75+shoulderOffs*p82+p78*p75+p54*p73+p60*p73-p75*p133-p78*p82+p82*p133-shoulderOffs*p76-shoulderOffs*p83+p66*p82+p42*p171+s3*lowArmLen*p30+s3*lowArmLen*p36+p24*p73+p56*p98+p28*p73+p77*p100+p77*p108-p49*p171-p61*p73+p78*p76+p78*p83+p63*p73-p76*p133-p83*p133-p66*p75-p51*p98-p24*p172-p56*p154-p28*p172-p77*p114-p66*p76-p66*p83-s3*lowArmLen*p31-p46*p98-p25*p73-p52*p98-p77*p115+p41*p171+p43*p171+s3*lowArmLen*p37+p29*p73+p77*p123-p51*p181-p42*p66+p51*p154+p56*p181+p49*p66+p77*p109-p46*p181-p52*p181-p41*p66-p43*p66+p46*p154+p25*p172+p52*p154+p77*p110+p77*p116-p29*p172+p42*p162-p49*p162+p41*p162+p43*p162;
    A[1,6] = -shoulderOffs*p75+shoulderOffs*p82+p78*p75+p54*p73+p60*p73-p75*p133-p78*p82+p82*p133-shoulderOffs*p76-shoulderOffs*p83+p66*p82+p42*p171+s3*lowArmLen*p30+s3*lowArmLen*p36+p24*p73+p56*p98+p28*p73+p77*p100+p77*p108-p49*p171-p61*p73+p78*p76+p78*p83+p63*p73-p76*p133-p83*p133-p66*p75-p51*p98-p24*p172-p56*p154-p28*p172-p77*p114+p49*p147-p66*p76-p66*p83-s3*lowArmLen*p31-p46*p98-p25*p73-p52*p98-p77*p115+p41*p171+p43*p171+s3*lowArmLen*p37+p29*p73+p77*p123-p51*p181-p42*p66+p51*p154+p100*p90+p90*p108+p56*p181+p49*p66+p77*p109-p42*p147-p30*p185-p36*p185-p90*p115-p46*p181-p52*p181-p41*p66-p43*p66+p31*p185+p46*p154+p25*p172+p52*p154+p90*p123+p77*p110+p77*p116-p41*p147-p43*p147-p37*p185-p29*p172+p51*p178+s2*p51*p57+p42*p162+p90*p109-s3*p24*p57-p56*p178-s2*p57*p56-s3*p28*p57-p49*p162-p90*p114+p46*p178+s2*p46*p57+s3*p25*p57+p52*p178+s2*p52*p57+p41*p162+p43*p162+p90*p110+p90*p116-s3*p29*p57+p4*p42*p57-p4*p49*p57+p4*p41*p57+p4*p43*p57;
    A[1,7] = -p75+p82-p76-p83;
    A[1,9] = +p54+p60-p61+p63;
    A[1,10] = -c2*p75+c2*p82-c2*p76-c2*p83;
    A[1,11] = +p54+p60-p61+p63;
    A[1,12] = -s2*p75+s2*p82-s2*p76-s2*p83;
    A[1,13] = +s3*p54+s3*p60+p4*p82-s3*p61+s3*p63-p4*p75-p4*p76-p4*p83;
    A[1,14] = +s2*p75-s2*p82+s2*p76+s2*p83;
    A[1,15] = +p100+p108+p109-p114-p115+p123+p110+p116;
    A[1,16] = +s2*p51-s3*p24-s2*p56-s3*p28+s2*p46+s3*p25+s2*p52-s3*p29+p4*p42-p4*p49+p4*p41+p4*p43;
    A[1,17] = +p100+p108+p109-p114-p115+p123+p110+p116;
    A[1,18] = +s2*p42+s3*p30+s3*p36-s2*p49-s3*p31+s2*p41+s2*p43+s3*p37-p4*p51+p4*p56-p4*p46-p4*p52;
    A[1,19] = +s5*p100+s5*p108-s5*p114-s5*p115+s5*p123+s2*c5*p51+s5*p109-p24*p119-s2*c5*p56-p28*p119+s2*c5*p46+p25*p119+s2*c5*p52+s5*p110+s5*p116-p29*p119+p42*p112-p49*p112+p41*p112+p43*p112;
    A[1,20] = -s2*p42-s3*p30-s3*p36+s2*p49+s3*p31-s2*p41-s2*p43-s3*p37+p4*p51-p4*p56+p4*p46+p4*p52;
    A[1,21] = +c5*p100+c5*p108+p24*p45+p56*p64+p28*p45-c5*p115+c5*p123+c5*p109-p51*p64-c5*p114+c5*p110+c5*p116-p46*p64-p25*p45-p52*p64+p29*p45-p42*p58+p49*p58-p41*p58-p43*p58;
    A[1,22] = +p19*p49+p18*p100+p18*p108-p19*p42-p20*p30-p20*p36-p18*p115+p20*p31+p18*p123-p19*p41-p19*p43-p20*p37+p23*p51+s2*p11*p51+p18*p109-s3*p11*p24-p23*p56-s2*p11*p56-s3*p11*p28-p18*p114+p23*p46+s2*p11*p46+s3*p11*p25+p23*p52+s2*p11*p52+p18*p110+p18*p116-s3*p11*p29+p4*p11*p42-p4*p11*p49+p4*p11*p41+p4*p11*p43;
    A[1,23] = +c5*p100+c5*p108+p24*p45+p56*p64+p28*p45-c5*p115+c5*p123+c5*p109-p51*p64-c5*p114+c5*p110+c5*p116-p46*p64-p25*p45-p52*p64+p29*p45-p42*p58+p49*p58-p41*p58-p43*p58;
    A[1,24] = -p49*p72+p42*p72+p30*p101+p36*p101+p53*p100+p53*p108-p31*p101-p53*p115+p41*p72+p43*p72+p37*p101+p53*p123-p51*p80+p47*p51+p56*p80+p53*p109-p24*p95-p47*p56-p28*p95-p53*p114-p46*p80-p52*p80+p47*p46+p25*p95+p47*p52+p53*p110+p53*p116-p29*p95+p42*p86-p49*p86+p41*p86+p43*p86;
    A[1,25] = +p70*p100+p70*p108-p70*p115+p70*p123+p70*p109+p71*p100+p71*p108-p51*p169-p42*p191-c7*p20*p30-c7*p20*p36-p70*p114+p24*p179+p56*p169+p28*p179+p49*p191-p71*p115+c7*p20*p31+p70*p110+p70*p116+p71*p123-p46*p169-p25*p179-p41*p191-p52*p169-p43*p191-c7*p20*p37+p29*p179+c7*p23*p51+s2*p44*p51+p71*p109-p42*p184-s3*p24*p44-c7*p23*p56-s2*p44*p56-s3*p28*p44-p71*p114+p49*p184+p71*p110+p71*p116-p41*p184-p43*p184-s3*p29*p44+p4*p44*p42-p4*p44*p49+c7*p23*p46+s2*p46*p44+s3*p25*p44+c7*p23*p52+s2*p44*p52+p4*p41*p44+p4*p44*p43;
    A[1,26] = +p65*p100+p65*p108-p49*p188-p65*p115+p65*p123+p65*p109-p51*p161-p79*p100-p65*p114-p79*p108+p42*p188+p24*p170+s7*p20*p30+p56*p161+p28*p170+s7*p20*p36+p65*p110+p79*p115+p65*p116-p46*p161-p25*p170-s7*p20*p31-p52*p161-p79*p123+p41*p188+p43*p188+p29*p170+s7*p20*p37-p42*p180-s7*p23*p51-s2*p51*p50-p79*p109+s3*p24*p50+p49*p180+s7*p23*p56+s2*p50*p56+s3*p28*p50+p79*p114-p41*p180-s7*p23*p46-s2*p46*p50-s3*p25*p50-p43*p180-s7*p23*p52-s2*p52*p50-p79*p110-p79*p116+s3*p29*p50-p4*p42*p50+p4*p50*p49-p4*p41*p50-p4*p43*p50;
    A[1,27] = -p49*p72+p42*p72+p30*p101+p36*p101+p53*p100+p53*p108-p31*p101-p53*p115+p41*p72+p43*p72+p37*p101+p53*p123-p51*p80+p47*p51+p56*p80+p53*p109-p24*p95-p47*p56-p28*p95-p53*p114-p46*p80-p52*p80+p47*p46+p25*p95+p47*p52+p53*p110+p53*p116-p29*p95+p42*p86-p49*p86+p41*p86+p43*p86;

    A[2,2] = -p5*p78-p54*p141-p60*p141+p61*p141-p63*p141;
    A[2,3] = -p5*p78+p5*p133+p5*p66-p54*p148-p60*p148-p54*p141-p60*p141+p54*p87+p60*p87+p61*p148+p61*p141-p63*p148-p63*p141-p61*p87+p63*p87;
    A[2,4] = -p5*p78+p5*p133+p5*p66+p32*p98-p54*p148-p60*p148-p54*p141-p60*p141+p54*p87+p60*p87+p61*p148+p61*p141-p63*p148-p63*p141+p26*p66-elbowOffs*p157-elbowOffs*p166+elbowOffs*p174-p61*p87-elbowOffs*p183+p63*p87+p24*p87+p28*p87-p25*p87+p29*p87;
    A[2,5] = -p5*p78+p5*p133+p5*p66+p32*p98-p26*p171-p54*p148-p60*p148-p54*p141-p60*p141+p54*p87+p60*p87-p77*p189+c2*lowArmLen*p24+c2*lowArmLen*p28+p61*p148+p61*p141-p63*p148-p63*p141+p32*p181+p26*p66-elbowOffs*p157-elbowOffs*p166-p32*p154-p77*p124-p77*p132+elbowOffs*p174-p61*p87-c2*lowArmLen*p25-elbowOffs*p183+p63*p87+c2*lowArmLen*p29-p26*p162+p38*p157+p38*p166+lowArmLen*p21*p30+lowArmLen*p21*p36+p24*p87+p28*p87-p38*p174-lowArmLen*p21*p31-p25*p87+p38*p183+lowArmLen*p21*p37+p29*p87+p77*p139-p77*p150-p38*p194-p38*p195+p38*p196-p38*p197;
    A[2,6] = -p5*p78+p5*p133+p5*p66+p32*p98-p26*p171-p54*p148-p60*p148-p54*p141-p60*p141+p54*p87+p60*p87-p77*p189+p26*p147+c2*lowArmLen*p24+c2*lowArmLen*p28+p61*p148+p61*p141-p63*p148-p63*p141+p32*p181+p26*p66-elbowOffs*p157-elbowOffs*p166-p32*p154-p77*p124-p77*p132+elbowOffs*p174-p61*p87-c2*lowArmLen*p25-elbowOffs*p183+p63*p87+c2*lowArmLen*p29-p32*p178-s2*p32*p57-p26*p162-p90*p189+p38*p157+p38*p166+lowArmLen*p21*p30+lowArmLen*p21*p36+p24*p87+p28*p87-handOffs*p24*p59-handOffs*p28*p59-p38*p174-lowArmLen*p21*p31-p25*p87+handOffs*p25*p59+p38*p183+lowArmLen*p21*p37+p29*p87+p77*p139-handOffs*p29*p59-p77*p150-p4*p26*p57+p57*p157+p57*p166-handOffs*p30*p55-handOffs*p36*p55-p38*p194-p38*p195-p90*p124-p90*p132+handOffs*p31*p55+p57*p183+p38*p196+p90*p139-handOffs*p37*p55-p38*p197-p90*p150-p57*p194-p57*p195-p57*p174-p57*p197+p57*p196;
    A[2,10] = +c2*p5+s2*p54+s2*p60-s2*p61+s2*p63;
    A[2,12] = +s2*p5-c2*p54-c2*p60+c2*p61-c2*p63;
    A[2,13] = +p5*p4+p21*p54+p21*p60-p21*p61+p21*p63;
    A[2,14] = -s2*p5+c2*p54+c2*p60-c2*p61+c2*p63;
    A[2,15] = -p189-p124-p132+p139-p150;
    A[2,16] = -s2*p32+p157+p166-p4*p26-p174+p183-p194-p195+p196-p197;
    A[2,17] = -p189-p124-p132+p139-p150;
    A[2,18] = -s2*p26+c2*p24+c2*p28+p4*p32-c2*p25+c2*p29+p21*p30+p21*p36-p21*p31+p21*p37;
    A[2,19] = -s5*p189-s2*c5*p32-s5*p124-s5*p132-p26*p112+c5*p157+c5*p166-c5*p174+c5*p183+s5*p139-s5*p150-c5*p194-c5*p195+c5*p196-c5*p197;
    A[2,20] = +s2*p26-c2*p24-c2*p28-p4*p32+c2*p25-c2*p29-p21*p30-p21*p36+p21*p31-p21*p37;
    A[2,21] = +p32*p64-c5*p189+p26*p58-s5*p157-s5*p166-c5*p124-c5*p132+s5*p174+c5*p139-s5*p183-c5*p150+s5*p194+s5*p195-s5*p196+s5*p197;
    A[2,22] = +p19*p26-p23*p32-s2*p11*p32-p18*p189-p24*p59-p28*p59+p25*p59-p29*p59-p4*p11*p26+p11*p157+p11*p166-p30*p55-p36*p55-p18*p124-p18*p132-p11*p174+p31*p55+p11*p183+p18*p139-p37*p55-p18*p150-p11*p194-p11*p195-p11*p197+p11*p196;
    A[2,23] = +p32*p64-c5*p189+p26*p58-s5*p157-s5*p166-c5*p124-c5*p132+s5*p174+c5*p139-s5*p183-c5*p150+s5*p194+s5*p195-s5*p196+s5*p197;
    A[2,24] = -p26*p72+p32*p80-p32*p47-p53*p189+c2*c6*p24+c2*c6*p28-c2*c6*p25+c2*c6*p29-p26*p86+c6*p21*p30+p30*p131+c6*p21*p36+p36*p131-p54*p186-p60*p186-c6*p21*p31-p31*p131+c6*p21*p37+p37*p131+p61*p186-p63*p186-p24*p173-p28*p173+p25*p173-p29*p173;
    A[2,25] = -p70*p189+p32*p169+p26*p191-c7*p23*p32-s2*p32*p44-p71*p189+p26*p184-p96*p157-c7*p24*p59-p96*p166-c7*p28*p59-p70*p124-p70*p132+p96*p174+c7*p25*p59+p70*p139-p96*p183-c7*p29*p59-p70*p150-p4*p26*p44+p44*p157+p44*p166-c7*p30*p55-c7*p36*p55-p71*p124-p71*p132+p96*p194+p96*p195+p71*p139-p96*p196-c7*p37*p55-p71*p150+p96*p197-p44*p194-p44*p195-p44*p174+c7*p31*p55+p44*p183+p44*p196-p44*p197;
    A[2,26] = -p65*p189+p32*p161-p26*p188+p26*p180+s7*p23*p32+s2*p32*p50+p79*p189-p88*p157-p88*p166-p65*p124-p65*p132+s7*p24*p59+s7*p28*p59+p88*p174+p65*p139-s7*p25*p59-p88*p183-p65*p150+s7*p29*p59+p4*p26*p50-p50*p157-p50*p166+p88*p194+s7*p30*p55+p88*p195+s7*p36*p55+p79*p124+p79*p132-p88*p196-s7*p31*p55-p50*p183-p79*p139+p88*p197+s7*p37*p55+p79*p150+p50*p194+p50*p195+p50*p174-p50*p196+p50*p197;
    A[2,27] = -p26*p72+p32*p80-p32*p47-p53*p189+c2*c6*p24+c2*c6*p28-c2*c6*p25+c2*c6*p29-p26*p86+c6*p21*p30+p30*p131+c6*p21*p36+p36*p131-p54*p186-p60*p186-c6*p21*p31-p31*p131+c6*p21*p37+p37*p131+p61*p186-p63*p186-p24*p173-p28*p173+p25*p173-p29*p173;

    A[3,3] = -elbowOffs*p48-elbowOffs*p67+elbowOffs*p74+elbowOffs*p84+elbowOffs*p93-elbowOffs*p62-elbowOffs*p68-elbowOffs*p85+elbowOffs*p94;
    A[3,4] = -elbowOffs*p67-elbowOffs*p48+elbowOffs*p74-elbowOffs*p117+elbowOffs*p84+elbowOffs*p93+elbowOffs*p145-elbowOffs*p62-elbowOffs*p68-elbowOffs*p137-elbowOffs*p85-elbowOffs*p155-elbowOffs*p163+elbowOffs*p94+elbowOffs*p130+elbowOffs*p146-elbowOffs*p156+elbowOffs*p164;
    A[3,5] = -elbowOffs*p67-elbowOffs*p48+elbowOffs*p74-lowArmLen*p32*p40-elbowOffs*p117+elbowOffs*p84+elbowOffs*p93-p77*p192+c3*lowArmLen*p56+elbowOffs*p145+p75*p89-p89*p82-elbowOffs*p62-elbowOffs*p68-c3*lowArmLen*p51-elbowOffs*p137-p38*p145-elbowOffs*p85-c3*lowArmLen*p46-c3*lowArmLen*p52-elbowOffs*p155-elbowOffs*p163+elbowOffs*p94+p76*p89+p89*p83+p38*p137+lowArmLen*p15*p30+lowArmLen*p15*p36+elbowOffs*p130+elbowOffs*p146+p38*p117+p54*p99+p60*p99-lowArmLen*p15*p31-elbowOffs*p156-p61*p99+p38*p155+p38*p163+lowArmLen*p15*p37+elbowOffs*p164+p63*p99-p38*p130-p38*p146+p38*p156-p38*p164;
    A[3,6] = -elbowOffs*p67-elbowOffs*p48+elbowOffs*p74-lowArmLen*p32*p40-elbowOffs*p117+elbowOffs*p84+elbowOffs*p93-p77*p192+c3*lowArmLen*p56+elbowOffs*p145+p75*p89-p89*p82-elbowOffs*p62-elbowOffs*p68-c3*lowArmLen*p51-elbowOffs*p137-p38*p145-s3*p82*p90-elbowOffs*p85-c3*lowArmLen*p46-c3*lowArmLen*p52-elbowOffs*p155-elbowOffs*p163+elbowOffs*p94+p76*p89+p89*p83-p90*p192+handOffs*p51*p182+p38*p137+s3*p32*p147+s3*p75*p90+lowArmLen*p15*p30+lowArmLen*p15*p36+elbowOffs*p130+elbowOffs*p146+p38*p117+p54*p99+p60*p99-handOffs*p56*p182-lowArmLen*p15*p31-elbowOffs*p156-p61*p99+handOffs*p46*p182+handOffs*p52*p182+p38*p155+p38*p163+s3*p76*p90+s3*p83*p90+lowArmLen*p15*p37+elbowOffs*p164+p63*p99+p57*p137+p57*p117+p4*p54*p90+p4*p60*p90-handOffs*p30*p187-handOffs*p36*p187-p57*p145-p38*p130-p38*p146-p4*p61*p90+handOffs*p31*p187+p57*p155+p57*p163+p38*p156+p4*p63*p90-handOffs*p37*p187-p38*p164-p57*p130-p57*p146-p57*p164+p57*p156;
    A[3,13] = -p67-p48+p74+p84+p93-p62-p68-p85+p94;
    A[3,15] = +s3*p75-p192-s3*p82+p4*p54+p4*p60+s3*p76+s3*p83-p4*p61+p4*p63;
    A[3,16] = +p137+p117-p145+p155+p163-p130-p146+p156-p164;
    A[3,17] = -p192+s3*p75-s3*p82+p4*p54+p4*p60+s3*p76+s3*p83-p4*p61+p4*p63;
    A[3,18] = -p32*p40+c3*p56-c3*p51-c3*p46-c3*p52+p15*p30+p15*p36-p15*p31+p15*p37;
    A[3,19] = -s5*p192+p45*p75-p45*p82+p45*p76+p45*p83+c5*p137+c5*p117+p54*p58+p60*p58-c5*p145-p61*p58+c5*p155+c5*p163+p58*p63-c5*p130-c5*p146+c5*p156-c5*p164;
    A[3,20] = +p32*p40-c3*p56+c3*p51+c3*p46+c3*p52-p15*p30-p15*p36+p15*p31-p15*p37;
    A[3,21] = +p75*p119-p82*p119-c5*p192+s5*p145+p76*p119+p83*p119+p54*p112+p60*p112-s5*p137-s5*p117-p61*p112+p63*p112-s5*p155-s5*p163+s5*p130+s5*p146-s5*p156+s5*p164;
    A[3,22] = -s3*p18*p82-p18*p192+p51*p182+s3*p19*p32+s3*p18*p75-p56*p182+p46*p182+p52*p182+s3*p18*p76+s3*p18*p83+p11*p137+p11*p117+p4*p18*p54+p4*p18*p60-p30*p187-p36*p187-p11*p145-p4*p18*p61+p31*p187+p11*p155+p11*p163+p4*p18*p63-p37*p187-p11*p130-p11*p146-p11*p164+p11*p156;
    A[3,23] = +p75*p119-p82*p119-c5*p192+s5*p145+p76*p119+p83*p119+p54*p112+p60*p112-s5*p137-s5*p117-p61*p112+p63*p112-s5*p155-s5*p163+s5*p130+s5*p146-s5*p156+s5*p164;
    A[3,24] = -p82*p138-c3*c6*p51-c6*p32*p40-s5*p5*p55+c3*c6*p56+p75*p138-c3*c6*p46-c3*c6*p52+p76*p138+p83*p138+c6*p15*p30+p102*p137+c6*p15*p36+s3*p26*p47+p54*p165+p60*p165-p102*p145-c6*p15*p31-p61*p165+p102*p155+c6*p15*p37+p102*p163+p63*p165-p102*p130-p102*p146+p102*p156-p102*p164;
    A[3,25] = -s3*p70*p82-p70*p192+s3*p75*p70+s3*p76*p70+s3*p70*p83-p71*p192+c7*p51*p182+s3*p32*p191+p4*p54*p70+s3*p75*p71+p4*p60*p70-p96*p137-c7*p56*p182-p96*p117-s3*p71*p82+p96*p145-p4*p61*p70+c7*p46*p182+c7*p52*p182+s3*p76*p71+p4*p63*p70+s3*p71*p83-p96*p155-p96*p163+p44*p137+p44*p117+p4*p54*p71+p4*p60*p71-c7*p30*p187-c7*p36*p187-p44*p145+p96*p130+p96*p146+p4*p63*p71-p96*p156-c7*p37*p187+p96*p164-p44*p130-p44*p146-p4*p61*p71+c7*p31*p187+p44*p155+p44*p163+p44*p156-p44*p164;
    A[3,26] = -p65*p192+s3*p65*p75-s3*p65*p82+s3*p65*p76+s3*p65*p83+p4*p54*p65+p4*p60*p65+p79*p192-p88*p137-s7*p51*p182-p88*p117-s3*p32*p188-s3*p75*p79+p88*p145+s7*p56*p182+s3*p82*p79-p4*p61*p65+p4*p65*p63-p88*p155-s7*p46*p182-p88*p163-s7*p52*p182-s3*p76*p79-s3*p83*p79-p50*p137-p50*p117-p4*p54*p79-p4*p60*p79+p88*p130+s7*p30*p187+p88*p146+s7*p36*p187+p50*p145-p88*p156-s7*p31*p187-p50*p155-p50*p163-p4*p63*p79+p88*p164+s7*p37*p187+p50*p130+p50*p146+p4*p61*p79-p50*p156+p50*p164;
    A[3,27] = -p82*p138-c3*c6*p51-c6*p32*p40-s5*p5*p55+c3*c6*p56+p75*p138-c3*c6*p46-c3*c6*p52+p76*p138+p83*p138+c6*p15*p30+p102*p137+c6*p15*p36+s3*p26*p47+p54*p165+p60*p165-p102*p145-c6*p15*p31-p61*p165+p102*p155+c6*p15*p37+p102*p163+p63*p165-p102*p130-p102*p146+p102*p156-p102*p164;

    A[4,4] = +elbowOffs*p104+p24*p98+p28*p98+p56*p73+p32*p87-p51*p73-p25*p98-p46*p73-p52*p73+p29*p98-p30*p66-p36*p66+p31*p66-p37*p66;
    A[4,5] = +lowArmLen*p69+elbowOffs*p104-lowArmLen*p97+lowArmLen*p106+lowArmLen*p127+lowArmLen*p113+p24*p98+p28*p98+p56*p73-lowArmLen*p134+p32*p87-p51*p73-p38*p104-p24*p154-p28*p154-p56*p172-lowArmLen*p120-p25*p98-p46*p73-p52*p73+lowArmLen*p142+lowArmLen*p125+lowArmLen*p152+p29*p98-p30*p66-p36*p66-p38*p107+p51*p172+lowArmLen*p121+lowArmLen*p135+p31*p66-lowArmLen*p126-p37*p66+p25*p154+p46*p172+p52*p172+lowArmLen*p143-p29*p154+p30*p162+p36*p162-p31*p162+p37*p162;
    A[4,6] = +lowArmLen*p69+elbowOffs*p104-lowArmLen*p97+lowArmLen*p106+lowArmLen*p127+lowArmLen*p113+p24*p98+p28*p98+p56*p73-lowArmLen*p134+p32*p87-handOffs*p32*p59-p51*p73-p38*p104-p24*p154-p28*p154-p56*p172+p49*p185-lowArmLen*p120-p25*p98-p46*p73-p52*p73+lowArmLen*p142+lowArmLen*p125+lowArmLen*p152+p29*p98-p57*p104-p30*p66-p36*p66-p38*p107+p51*p172+handOffs*p151+lowArmLen*p121+lowArmLen*p135-p30*p147-p42*p185-p36*p147+p31*p66-lowArmLen*p126-p37*p66+p31*p147+p25*p154+p46*p172+p52*p172+lowArmLen*p143-p41*p185-p37*p147-p43*p185-p29*p154-p57*p107+s3*p51*p57+p30*p162+p36*p162-p24*p178-s2*p24*p57-p28*p178-s2*p28*p57-s3*p57*p56-p31*p162+p25*p178+s2*p25*p57+s3*p46*p57+s3*p52*p57+p37*p162-p29*p178-s2*p29*p57+p4*p30*p57+p4*p36*p57-p4*p31*p57+p4*p37*p57;
    A[4,16] = -p104+s3*p51-s2*p24-s2*p28-s3*p56-p107+s2*p25+s3*p46+s3*p52-s2*p29+p4*p30+p4*p36-p4*p31+p4*p37;
    A[4,18] = +p69-p97+p106+p127+p113-p134-p120+p142+p125+p152+p121+p135-p126+p143;
    A[4,19] = -c5*p104-c5*p107+p51*p119-s2*c5*p24-s2*c5*p28-p56*p119+s2*c5*p25+p46*p119+p52*p119-s2*c5*p29+p30*p112+p36*p112-p31*p112+p37*p112;
    A[4,20] = -p69+p97-p106-p127-p113+p134+p120-p142-p125-p152-p121-p135+p126-p143;
    A[4,21] = +s5*p104+p24*p64+p28*p64+p45*p56+s5*p107-p45*p51-p25*p64-p46*p45-p45*p52+p29*p64-p30*p58-p36*p58+p31*p58-p37*p58;
    A[4,22] = -p32*p59+p20*p49-p11*p104+p151-p19*p30-p20*p42-p19*p36+p19*p31-p20*p41-p19*p37-p20*p43-p11*p107+s3*p11*p51-p24*p23-s2*p11*p24-p23*p28-s2*p11*p28-s3*p11*p56+p25*p23+s2*p11*p25+s3*p11*p46+s3*p11*p52-p23*p29-s2*p11*p29+p4*p11*p30+p4*p11*p36-p4*p11*p31+p4*p11*p37;
    A[4,23] = +s5*p104+p24*p64+p28*p64+p45*p56+s5*p107-p45*p51-p25*p64-p46*p45-p45*p52+p29*p64-p30*p58-p36*p58+p31*p58-p37*p58;
    A[4,24] = +c6*p69-p49*p101-c6*p97-p26*p131+p30*p72+p42*p101+p36*p72-p31*p72+p41*p101+p37*p72+p43*p101-p32*p173+p24*p80+p51*p95+p28*p80-p24*p47-p28*p47-p56*p95-p25*p80+p25*p47+p46*p95+p29*p80+p52*p95-p29*p47+p30*p86+p36*p86+p37*p86-p31*p86;
    A[4,25] = -c7*p32*p59+p96*p104-p44*p104+p96*p107+c7*p151-p51*p179-p30*p191-c7*p20*p42-p36*p191+p24*p169+p28*p169+p56*p179+c7*p20*p49+p31*p191-p25*p169-p46*p179-c7*p20*p41-p52*p179-p37*p191-c7*p20*p43+p29*p169-p44*p107+s3*p44*p51-p30*p184-c7*p24*p23-s2*p24*p44-p36*p184-c7*p23*p28-s2*p28*p44-s3*p44*p56-p37*p184-c7*p23*p29-s2*p29*p44+p4*p30*p44+p4*p36*p44+p31*p184+c7*p25*p23+s2*p25*p44+s3*p46*p44+s3*p44*p52-p4*p31*p44+p4*p37*p44;
    A[4,26] = +p88*p104+s7*p32*p59-s7*p20*p49+p88*p107+p50*p104-p51*p170-s7*p151+p24*p161+p30*p188+s7*p20*p42+p28*p161+p36*p188+p56*p170-p25*p161-p31*p188-p46*p170-p52*p170+s7*p20*p41+p29*p161+p37*p188+s7*p20*p43+p50*p107-p30*p180-s3*p51*p50-p36*p180+s7*p24*p23+s2*p24*p50+s7*p23*p28+s2*p28*p50+s3*p50*p56-s7*p25*p23-s2*p25*p50-s3*p46*p50-p37*p180-s3*p52*p50+s7*p23*p29+s2*p29*p50-p4*p30*p50-p4*p36*p50+p31*p180+p4*p31*p50-p4*p37*p50;
    A[4,27] = +c6*p69-p49*p101-c6*p97-p26*p131+p30*p72+p42*p101+p36*p72-p31*p72+p41*p101+p37*p72+p43*p101-p32*p173+p24*p80+p51*p95+p28*p80-p24*p47-p28*p47-p56*p95-p25*p80+p25*p47+p46*p95+p29*p80+p52*p95-p29*p47+p30*p86+p36*p86+p37*p86-p31*p86;

    A[5,5] = -p48*p38+p38*p74-p69*p77-p38*p67-wristOffs*p30*p64-p42*p89-wristOffs*p36*p64+p49*p89-p38*p62-p38*p68+p77*p97+p38*p84+p38*p93-p38*p85+wristOffs*p31*p64+p38*p94-p41*p89-wristOffs*p37*p64-p43*p89-p24*p99-p28*p99+p25*p99-p29*p99;
    A[5,6] = -p48*p38+p38*p74-p69*p77-p38*p67-wristOffs*p30*p64-p42*p89-wristOffs*p36*p64+p49*p89-p38*p62-p38*p68-p69*p90-p57*p67-p48*p57+p57*p74+p77*p97+p38*p84+p38*p93-p62*p57-p57*p68-p38*p85+wristOffs*p31*p64+p38*p94-p41*p89-wristOffs*p37*p64-p43*p89+p97*p90+p57*p84+p57*p93-p90*p106-p90*p127-p90*p113-p24*p99-p28*p99+p90*p134-p57*p85+p90*p120+p25*p99+p57*p94-p90*p142-p90*p125-p90*p152-p29*p99-p90*p121-p90*p135-p90*p143+p90*p126;
    A[5,19] = -c5*p48+c5*p74-s5*p69-c5*p67+p45*p49-c5*p62-c5*p68+s5*p97+c5*p84+c5*p93-p30*p64-p45*p42-p36*p64-c5*p85+p31*p64+c5*p94-p41*p45-p37*p64-p45*p43-p24*p58-p28*p58+p25*p58-p29*p58;
    A[5,21] = +s5*p67+s5*p48-s5*p74-c5*p69+p49*p119+s5*p62+s5*p68+c5*p97-c5*p106-p42*p119-c5*p113-s5*p84-s5*p93+c5*p120+s5*p85-p41*p119-c5*p125-p43*p119-s5*p94-p24*p112-p28*p112+p25*p112-p29*p112;
    A[5,22] = -p18*p69-p11*p67-p11*p48+p11*p74-p11*p62-p11*p68+p18*p97+p11*p84+p11*p93-p18*p106-p18*p127-p18*p113+p18*p134-p11*p85+p18*p120+p11*p94-p18*p142-p18*p125-p18*p152-p18*p121-p18*p135-p18*p143+p18*p126;
    A[5,23] = +s5*p67+s5*p48-s5*p74-c5*p69+p49*p119+s5*p62+s5*p68+c5*p97-c5*p106-p42*p119-c5*p113-s5*p84-s5*p93+c5*p120+s5*p85-p41*p119-c5*p125-p43*p119-s5*p94-p24*p112-p28*p112+p25*p112-p29*p112;
    A[5,24] = -p53*p69-p67*p102-s3*p5*p47+p74*p102+p49*p138-p62*p102-p68*p102+s5*p151+p84*p102+p93*p102-s5*p19*p30-p42*p138-s5*p19*p36-p85*p102+s5*p19*p31+p94*p102-p41*p138-s5*p19*p37-p43*p138-p24*p165-p28*p165+p25*p165-p29*p165;
    A[5,25] = -p74*p96-p69*p70+p67*p96+p48*p96+p62*p96+p68*p96-p69*p71-p44*p67-p48*p44+p70*p97+p44*p74-p70*p106-p70*p127-p70*p113-p84*p96-p96*p93+p70*p134-p44*p62-p44*p68+p70*p120+p85*p96-p70*p142-p70*p125-p70*p152-p96*p94+p71*p97+p44*p84+p44*p93-p70*p121-p71*p106-p71*p127-p70*p135-p71*p113+p71*p134+p44*p94-p71*p142-p70*p143-p71*p125-p71*p152-p71*p121-p71*p135-p44*p85+p70*p126+p71*p120+p71*p126-p71*p143;
    A[5,26] = -p69*p65+p67*p88+p48*p88-p74*p88+p62*p88+p68*p88+p65*p97+p69*p79+p50*p67+p48*p50-p65*p106-p65*p127-p65*p113-p88*p84-p88*p93-p50*p74+p65*p134+p65*p120+p88*p85+p50*p62+p50*p68-p65*p142-p65*p125-p65*p152-p88*p94-p65*p121-p65*p135-p79*p97-p50*p84-p50*p93+p79*p106+p79*p127+p79*p113-p79*p134-p79*p120-p65*p143-p50*p94+p79*p142+p79*p125+p79*p152+p79*p121+p79*p135+p65*p126+p50*p85-p79*p126+p79*p143;
    A[5,27] = -p53*p69-p67*p102-s3*p5*p47+p74*p102+p49*p138-p62*p102-p68*p102+s5*p151+p84*p102+p93*p102-s5*p19*p30-p42*p138-s5*p19*p36-p85*p102+s5*p19*p31+p94*p102-p41*p138-s5*p19*p37-p43*p138-p24*p165-p28*p165+p25*p165-p29*p165;

    A[6,6] = -c6*p104*p158-c6*p107*p158-p32*p105*p131-c6*p105*p104+p51*p101*p158+p53*p67*p105+p5*p20*p53*p147-p24*p72*p158-p28*p72*p158-p56*p101*p158-p53*p74*p105-handOffs*p91*p118+p30*p80*p158+p36*p80*p158+p26*p105*p173+p51*p91*p185+p67*p91*p90+p5*p22*p45*p147-p30*p47*p105-p42*p95*p105-p36*p47*p105-p53*p84*p105-p24*p91*p147-p53*p93*p105-p28*p91*p147-p56*p91*p185-p74*p91*p90+p49*p95*p105+p25*p72*p158+p46*p101*p158+p52*p101*p158+p53*p62*p105+p53*p68*p105-p29*p72*p158-p32*p59*p159-p84*p91*p90-p91*p93*p90+p20*p49*p159-p31*p80*p158+p31*p47*p105+p37*p80*p158+p53*p85*p105+p25*p91*p147+p46*p91*p185+p62*p91*p90+p52*p91*p185+p68*p91*p90-p41*p95*p105-p37*p47*p105-p43*p95*p105-p53*p94*p105-p29*p91*p147+p151*p159+p30*p91*p178+p36*p91*p178-p24*p86*p105-p19*p30*p159-p20*p42*p159-p28*p86*p105-p19*p36*p159-p31*p91*p178+p25*p86*p105+p19*p31*p159+p85*p91*p90+p37*p91*p178-p20*p41*p159-p29*p86*p105-p19*p37*p159-p20*p43*p159-p91*p94*p90-p23*p24*p159-p23*p28*p159+p23*p25*p159-p23*p29*p159;
    A[6,22] = +c6*p22*p104-c6*p104-c6*p104-p53*p74-p53*p74+c6*p22*p107+p32*p35*p131+c6*p35*p104-p22*p51*p101-p35*p53*p67-p5*p19*p20*p53+p22*p24*p72+p22*p28*p72+p22*p56*p101+p35*p53*p74-c6*p107-c6*p107-p32*p131-p32*p131+p51*p101+p51*p101+p53*p67+p53*p67+p5*p186+p5*p186-p24*p72-p24*p72-p28*p72-p56*p101-p28*p72-p56*p101+p32*p59*p149+p91*p118-p22*p30*p80-p22*p36*p80-p26*p35*p173-p20*p51*p91-p18*p67*p91-p5*p22*p186+p30*p35*p47+p35*p42*p95+p35*p36*p47+p35*p53*p84+p19*p24*p91+p35*p53*p93+p19*p28*p91+p20*p56*p91+p18*p74*p91+p25*p72+p46*p101+p25*p72+p46*p101+p52*p101+p52*p101-p35*p49*p95+p53*p62+p53*p62+p53*p68+p53*p68-p29*p72-p29*p72-p22*p25*p72-p22*p46*p101-p22*p52*p101+p30*p80+p30*p80+p36*p80+p36*p80-p35*p53*p62-p35*p53*p68+p26*p173+p26*p173+p22*p29*p72-p30*p47-p42*p95-p30*p47-p42*p95-p36*p47-p36*p47-p53*p84-p53*p84-p53*p93-p53*p93+p49*p95+p49*p95-p31*p80-p31*p80+p24*p35*p86+p19*p30*p149+p20*p42*p149+p28*p35*p86+p19*p36*p149+p18*p84*p91+p18*p93*p91+p31*p47+p31*p47+p37*p80+p37*p80-p20*p49*p149+p53*p85+p53*p85-p41*p95-p41*p95-p37*p47-p43*p95-p37*p47-p43*p95-p53*p94-p53*p94+p22*p31*p80-p31*p35*p47-p22*p37*p80-p35*p53*p85-p19*p25*p91-p20*p46*p91-p18*p62*p91-p20*p52*p91-p18*p68*p91+p35*p41*p95+p35*p37*p47+p35*p43*p95-p24*p86-p24*p86-p28*p86-p28*p86+p35*p53*p94+p19*p29*p91-p151*p149-p23*p30*p91-p23*p36*p91-p29*p86-p29*p86+p23*p31*p91-p25*p35*p86-p19*p31*p149-p18*p85*p91-p23*p37*p91+p20*p41*p149+p29*p35*p86+p19*p37*p149+p20*p43*p149+p18*p94*p91+p23*p24*p149+p23*p28*p149+p25*p86+p25*p86-p23*p25*p149+p23*p29*p149;
    A[6,24] = -p140-p140-p19*p24-p19*p24-p19*p28-p20*p56-p19*p28-p20*p56+p16*p35*p140+p17*p35*p140+p11*p69+p11*p69-p118-p118-p18*p67-p18*p67-p18*p48-p18*p48+p20*p51+p20*p51+p18*p74+p18*p74-p69*p167-p18*p62-p18*p62-p18*p68-p18*p68-p69*p175+p19*p25+p20*p46+p19*p25+p20*p46+p20*p52+p20*p52+p16*p35*p118+p16*p22*p140+p22*p67*p103+p17*p35*p118+p22*p48*p103+p17*p22*p140-p35*p51*p122+p22*p67*p111-p19*p29-p19*p29+p22*p48*p111-p35*p51*p128-p22*p74*p103+p24*p35*p129+p28*p35*p129+p35*p56*p122-p22*p74*p111+p24*p35*p136+p28*p35*p136+p35*p56*p128-p11*p97-p11*p97+p23*p30+p23*p30+p11*p106+p11*p127+p11*p106+p11*p127+p23*p36+p23*p36+p11*p113+p11*p113+p18*p84+p18*p84+p18*p93+p18*p93-p11*p134-p11*p134-p25*p35*p136-p35*p46*p128-p35*p52*p128+p29*p35*p129-p23*p31-p23*p31-p11*p120-p11*p120+p29*p35*p136-p18*p85-p18*p85-p35*p69*p92+p11*p142+p11*p142+p23*p37+p23*p37+p11*p125+p11*p152+p11*p125+p11*p152+p16*p22*p118+p97*p167-p35*p69*p81+p18*p94+p18*p94+p17*p22*p118+p97*p175-p30*p35*p153-p106*p167-p22*p51*p122-p127*p167+p35*p67*p103-p35*p36*p153-p113*p167+p5*p19*p20*p103-p30*p35*p144-p22*p84*p103-p106*p175-p22*p51*p128-p127*p175+p35*p67*p111-p35*p36*p144-p22*p93*p103-p113*p175+p5*p19*p20*p111-p22*p84*p111+p22*p24*p129-p22*p93*p111+p22*p28*p129+p22*p56*p122+p134*p167-p35*p74*p103+p22*p24*p136+p22*p28*p136+p22*p56*p128+p134*p175-p35*p74*p111+p11*p121+p11*p121+p11*p135+p11*p135+p22*p62*p103+p22*p68*p103-p25*p35*p129-p35*p46*p122+p22*p62*p111-p35*p52*p122+p22*p68*p111+p11*p143+p11*p143+p151*p190-p22*p30*p153-p121*p167-p22*p36*p153-p135*p167+p151*p193-p22*p30*p144-p121*p175-p19*p30*p190-p20*p42*p190-p22*p36*p144-p135*p175-p19*p36*p190-p20*p54*p59*p103-p19*p30*p193-p20*p42*p193-p20*p59*p60*p103-p19*p36*p193-p20*p54*p59*p111-p20*p59*p60*p111+p20*p49*p190+p20*p49*p193+p31*p35*p153+p120*p167+p31*p35*p144+p22*p85*p103+p120*p175+p22*p85*p111-p22*p25*p129-p22*p46*p122-p142*p167+p35*p62*p103-p35*p37*p153-p125*p167-p22*p52*p122-p152*p167+p35*p68*p103-p22*p25*p136-p22*p46*p128-p142*p175+p35*p62*p111-p35*p37*p144-p22*p94*p103-p125*p175-p22*p52*p128-p152*p175+p35*p68*p111-p22*p94*p111+p22*p29*p129-p11*p126-p11*p126+p22*p29*p136-p35*p92*p121-p23*p28*p190-p35*p81*p121-p23*p28*p193+p22*p31*p153+p126*p167+p22*p31*p144+p126*p175+p19*p31*p190+p20*p59*p61*p103+p19*p31*p193-p22*p37*p153-p143*p167+p20*p59*p61*p111-p20*p41*p190-p22*p37*p144-p143*p175-p19*p37*p190-p20*p43*p190-p20*p41*p193-p20*p59*p63*p103-p19*p37*p193-p20*p43*p193-p20*p59*p63*p111+p35*p92*p126+p35*p81*p126-p23*p29*p190-p23*p29*p193;
    A[6,25] = -c6*c7*p104-c6*c7*p107-p32*p59*p65+c7*p51*p101+c7*p53*p67+p5*p19*p170-c7*p24*p72-c7*p28*p72-c7*p56*p101-c7*p53*p74+c7*p25*p72+c7*p46*p101+c7*p52*p101+c7*p53*p62+c7*p53*p68-c7*p29*p72+c7*p30*p80+c7*p36*p80+c7*p26*p173-c7*p30*p47-p20*p42*p65-c7*p36*p47-c7*p53*p84-c7*p53*p93+p20*p49*p65+c7*p53*p85-p20*p41*p65-c7*p37*p47-p20*p43*p65-c7*p53*p94-p23*p24*p65-p23*p28*p65-c7*p31*p80+c7*p31*p47+c7*p37*p80+p23*p25*p65-p23*p29*p65;
    A[6,26] = +s7*c6*p104+s7*c6*p107+p32*p59*p70-s7*p51*p101-s7*p53*p67-p5*p19*p179+s7*p24*p72+s7*p28*p72+s7*p56*p101+s7*p53*p74-s7*p25*p72-s7*p46*p101-s7*p52*p101-s7*p53*p62-s7*p53*p68+s7*p29*p72-s7*p30*p80-s7*p36*p80-s7*p26*p173+s7*p30*p47+p20*p42*p70+s7*p36*p47+s7*p53*p84+s7*p53*p93-p20*p49*p70-s7*p31*p47-s7*p37*p80-s7*p53*p85+p20*p41*p70+s7*p37*p47+p20*p43*p70+s7*p53*p94+p23*p24*p70+p23*p28*p70+s7*p31*p80-p23*p25*p70+p23*p29*p70;
    A[6,27] = -p16*p140-p17*p140+p69*p81-p16*p118-p67*p103-p17*p118-p48*p103+p51*p122-p67*p111-p48*p111+p51*p128+p74*p103-p24*p129-p28*p129-p56*p122+p74*p111-p24*p136-p28*p136-p56*p128+p69*p92+p30*p144+p84*p103+p81*p106+p81*p127+p36*p144+p93*p103+p81*p113+p84*p111+p93*p111-p92*p134-p81*p134-p62*p103-p68*p103+p25*p129+p46*p122-p62*p111+p52*p122-p68*p111+p25*p136+p46*p128+p52*p128-p29*p129-p29*p136-p97*p92-p81*p97+p30*p153+p92*p106+p92*p127+p36*p153+p92*p113-p31*p153-p92*p120-p31*p144-p85*p103-p81*p120-p85*p111+p92*p142+p37*p153+p92*p125+p92*p152+p81*p142+p37*p144+p94*p103+p81*p125+p81*p152+p94*p111+p92*p121+p92*p135+p81*p121+p81*p135-p92*p126-p81*p126+p92*p143+p81*p143;

    A[7,25] = +s7*p17*p140-p69*p88-p67*p65-p48*p65+p65*p74-p69*p160+s7*p17*p118+s7*p16*p140+p17*p67*p79+p17*p48*p79-s7*p51*p128-p62*p65-p68*p65-p17*p74*p79+s7*p24*p136+s7*p28*p136+s7*p56*p128+p88*p97+p65*p84+p65*p93-p30*p161-p42*p170-p36*p161+p49*p170+p16*p48*p79-s7*p30*p144-p106*p160-p127*p160-s7*p51*p122-s7*p36*p144-p113*p160-p65*p85-p17*p84*p79-p17*p79*p93-p16*p74*p79+p31*p161+s7*p24*p129+s7*p28*p129+p134*p160+s7*p56*p122+p65*p94-p41*p170-p37*p161-p43*p170+p17*p62*p79+p17*p68*p79-s7*p25*p136-s7*p46*p128-s7*p52*p128-p24*p180-p28*p180+s7*p29*p136-p69*p176+p97*p160+s7*p16*p118+p16*p67*p79-p29*p180+s7*p31*p144+p120*p160+p17*p85*p79+p16*p62*p79+p16*p68*p79-s7*p25*p129-p142*p160-s7*p46*p122-s7*p37*p144-p125*p160-p160*p152-s7*p52*p122-p17*p79*p94+s7*p29*p129+p97*p176-p121*p160-s7*p30*p153-p106*p176-p127*p176-p135*p160-s7*p36*p153-p113*p176-p16*p84*p79-p16*p79*p93+p25*p180+p134*p176+p126*p160+s7*p31*p153+p120*p176+p16*p85*p79-p142*p176-p143*p160-s7*p37*p153-p125*p176-p152*p176-p16*p79*p94-p121*p176-p135*p176+p126*p176-p143*p176;
    A[7,26] = +c7*p16*p140+p69*p96+p67*p70+p48*p70-p74*p70-p49*p179+c7*p16*p118+p16*p67*p71+p16*p48*p71+c7*p17*p140-c7*p51*p122-p16*p74*p71+c7*p24*p129+c7*p28*p129+c7*p56*p122+p62*p70+p68*p70-p97*p96-p70*p84-p70*p93+p30*p169+p42*p179+p36*p169-p69*p168+p17*p48*p71-p16*p71*p84-c7*p51*p128-p16*p71*p93+p134*p168-p17*p74*p71+p70*p85+c7*p24*p136+c7*p28*p136+c7*p56*p128-p31*p169-p70*p94+p41*p179+p37*p169+p43*p179+p16*p62*p71+p16*p68*p71-c7*p25*p129-c7*p46*p122-c7*p52*p122+c7*p29*p129+p24*p184+p28*p184-p69*p177+p97*p168+c7*p17*p118-c7*p30*p153-p106*p168-p127*p168+p17*p67*p71-c7*p36*p153-p113*p168-p25*p184+p29*p184+c7*p31*p153+p120*p168+p16*p71*p85-p142*p168+p17*p62*p71-c7*p37*p153-p125*p168-p152*p168+p17*p68*p71-c7*p25*p136-c7*p46*p128-p16*p71*p94-c7*p52*p128+c7*p29*p136+p97*p177-p121*p168-p135*p168-c7*p30*p144-p106*p177-p127*p177-c7*p36*p144-p113*p177-p17*p71*p84-p17*p71*p93+p134*p177+p126*p168+c7*p31*p144+p120*p177+p17*p71*p85-p143*p168-p142*p177-c7*p37*p144-p125*p177-p152*p177-p17*p71*p94-p121*p177-p135*p177+p126*p177-p143*p177;
    A = A[:,[1:7; 9:27]];
    return A
end


text:
@testset "choice_at combinator" begin

    at = choice_at(bernoulli, Int)

    @testset "assess" begin
        choices = choicemap()
        choices[3] = true
        (weight, value) = assess(at, (0.4, 3), choices)
        @test isapprox(weight, log(0.4))
        @test value == true
    end

    @testset "propose" begin
        (choices, weight, value) = propose(at, (0.4, 3))
        @test isapprox(weight, value ? log(0.4) : log(0.6))
        @test choices[3] == value
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
    end

    @testset "generate" begin

        # without assignment
        (trace, weight) = generate(at, (0.4, 3), EmptyChoiceMap())
        @test weight == 0.
        value = get_retval(trace)
        choices = get_choices(trace)
        @test choices[3] == value
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0

        # with constraints
        constraints = choicemap()
        constraints[3] = true
        (trace, weight) = generate(at, (0.4, 3), constraints)
        value = get_retval(trace)
        @test value == true
        @test isapprox(weight, log(0.4))
        choices = get_choices(trace)
        @test choices[3] == value
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
    end

    function get_trace()
        constraints = choicemap()
        constraints[3] = true
        (trace, _) = generate(at, (0.4, 3), constraints)
        trace
    end

    @testset "project" begin
        trace = get_trace()
        @test isapprox(project(trace, EmptySelection()), 0.)
        selection = select(3)
        @test isapprox(project(trace, selection), log(0.4))
    end

    @testset "update" begin
        trace = get_trace()

        # change kernel_args, same key, no constraint
        (new_trace, weight, retdiff, discard) = update(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), EmptyChoiceMap())
        choices = get_choices(new_trace)
        @test choices[3] == true
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test isapprox(weight, log(0.2) - log(0.4))
        @test get_retval(new_trace) == true
        @test isempty(discard)

        # change kernel_args, same key, with constraint
        constraints = choicemap()
        constraints[3] = false
        (new_trace, weight, retdiff, discard) = update(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(new_trace)
        @test choices[3] == false
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test isapprox(weight, log(1 - 0.2) - log(0.4))
        @test get_retval(new_trace) == false
        @test discard[3] == true
        @test length(collect(get_values_shallow(discard))) == 1
        @test length(collect(get_submaps_shallow(discard))) == 0

        # change kernel_args, different key, with constraint
        constraints = choicemap()
        constraints[4] = false
        (new_trace, weight, retdiff, discard) = update(trace,
            (0.2, 4), (UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(new_trace)
        @test choices[4] == false
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test isapprox(weight, log(1 - 0.2) - log(0.4))
        @test get_retval(new_trace) == false
        @test discard[3] == true
        @test length(collect(get_values_shallow(discard))) == 1
        @test length(collect(get_submaps_shallow(discard))) == 0
    end

    @testset "regenerate" begin
        trace = get_trace()

        # change kernel_args, same key, not selected
        (new_trace, weight, retdiff) = regenerate(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), EmptySelection())
        choices = get_choices(new_trace)
        @test choices[3] == true
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test isapprox(weight, log(0.2) - log(0.4))
        @test get_retval(new_trace) == true
        @test isapprox(get_score(new_trace), log(0.2))

        # change kernel_args, same key, selected
        selection = select(3)
        (new_trace, weight, retdiff) = regenerate(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), selection)
        choices = get_choices(new_trace)
        value = choices[3]
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test weight == 0.
        @test get_retval(new_trace) == value
        @test isapprox(get_score(new_trace), log(value ? 0.2 : 1 - 0.2))

        # change kernel_args, different key, not selected
        (new_trace, weight, retdiff) = regenerate(trace,
            (0.2, 4), (UnknownChange(), UnknownChange()), EmptySelection())
        choices = get_choices(new_trace)
        value = choices[4]
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test weight == 0.
        @test get_retval(new_trace) == value
        @test isapprox(get_score(new_trace), log(value ? 0.2 : 1 - 0.2))
    end

    @testset "choice_gradients" begin
        y = 1.2
        constraints = choicemap()
        set_value!(constraints, 3, y)
        (trace, _) = generate(choice_at(normal, Int), (0.0, 1.0, 3), constraints)

        # not selected
        (input_grads, choices, gradients) = choice_gradients(
            trace, EmptySelection(), 1.2)
        @test isempty(choices)
        @test isempty(gradients)
        @test length(input_grads) == 3
        @test isapprox(input_grads[1], logpdf_grad(normal, y, 0.0, 1.0)[2])
        @test isapprox(input_grads[2], logpdf_grad(normal, y, 0.0, 1.0)[3])
        @test input_grads[3] == nothing # the key has no gradient

        # selected with retval_grad
        retval_grad = 1.234
        selection = select(3)
        (input_grads, choices, gradients) = choice_gradients(
            trace, selection, retval_grad)
        @test choices[3] == y
        @test isapprox(gradients[3], logpdf_grad(normal, y, 0.0, 1.0)[1] + retval_grad)
        @test length(collect(get_values_shallow(gradients))) == 1
        @test length(collect(get_submaps_shallow(gradients))) == 0
        @test length(collect(get_values_shallow(choices))) == 1
        @test length(collect(get_submaps_shallow(choices))) == 0
        @test length(input_grads) == 3
        @test isapprox(input_grads[1], logpdf_grad(normal, y, 0.0, 1.0)[2])
        @test isapprox(input_grads[2], logpdf_grad(normal, y, 0.0, 1.0)[3])
        @test input_grads[3] == nothing # the key has no gradient
    end

    @testset "accumulate_param_gradients!" begin
        trace = get_trace()
        input_grads = accumulate_param_gradients!(trace, nothing)
        @test length(input_grads) == 2
        @test isapprox(input_grads[1], logpdf_grad(bernoulli, true, 0.4)[2])
        @test input_grads[2] == nothing # the key has no gradient
    end

end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

import Pluto, MLCourse

session = Pluto.ServerSession()
session.options.server.port = 40404
session.options.security.require_secret_for_access = false

# for file in last.(MLCourse.NOTEBOOKS)
#     Base.include(Module(), joinpath(@__DIR__, "..", "notebooks", file))
# end
using CSV, DataFrames, MLJ, MLJLinearModels, NearestNeighborModels,
      Random, Distributions, MLJGLMInterface, Plots, StatsPlots, Statistics,
      LinearAlgebra

weather = CSV.read(joinpath(@__DIR__, "..", "data", "weather2015-2018.csv"), DataFrame)
@df weather corrplot([:BAS_pressure :LUG_pressure :LUZ_pressure :LUZ_wind_peak],
                     grid = false, fillcolor = cgrad(), size = (700, 600))

X = float.(select(weather, Not(:LUZ_wind_peak)))[1:1000, :]
y = weather.LUZ_wind_peak[1:1000]

m = fit!(machine(MLJLinearModels.LinearRegressor(), X, y))
predict(m, X)

m = fit!(machine(MLJGLMInterface.LinearRegressor(), X, y))
predict(m, X)
predict_mean(m, X)

pdf(Normal(0., 1.), 2.)
pdf(Bernoulli(.3), 1)
pdf(Categorical([.2, .3, .5]), [1, 0, 0])

for backend in (gr(), plotly())
    scatter(rand(10), rand(10))
    plot(sin)
    contour(0:10, 0:10, (x, y) -> tanh(x + y))
end

spamdata = CSV.read(joinpath(@__DIR__, "..", "data", "spam.csv"), DataFrame)
dropmissing!(spamdata)
import TextAnalysis: Corpus, StringDocument, DocumentTermMatrix, lexicon,
                     update_lexicon!, tf
crps = Corpus(StringDocument.(spamdata.text[1:200]))
update_lexicon!(crps)
lexicon(crps)
small_lex = Dict(k => lexicon(crps)[k]
                 for k in findall(x -> 80 <= last(x) <= 10^2, lexicon(crps)))
m = DocumentTermMatrix(crps, small_lex)
spam_or_ham = coerce(String.(spamdata.label[1:200]), Binary)
normalized_word_counts = float.(DataFrame(tf(m), :auto))
m3 = fit!(machine(LinearBinaryClassifier(),
                  normalized_word_counts,
                  spam_or_ham));
confusion_matrix(predict_mode(m3, normalized_word_counts), spam_or_ham)
function losses(machine, input, response)
    (logikelihood = -sum(log_loss(predict(machine, input), response)),
     misclassification_rate = mean(predict_mode(machine, input) .!= response),
     accuracy = accuracy(predict_mode(machine, input), response),
     auc = auc(predict(machine, input), response))
end;
losses(m3, normalized_word_counts, spam_or_ham)

begin
    f(x) = .3 * sin(10x) + .7x
    # f(x) = .0015*(12x-6)^4 -.035(12x-6)^2 + .5x + .2  # an alternative
	σ(x) = 1 / (1 + exp(-x))
    function regression_data_generator(; n, rng = Random.GLOBAL_RNG)
        x = rand(rng, n)
        DataFrame(x = x, y = f.(x) .+ .1*randn(rng, n))
    end
    function classification_data_generator(; n, rng = Random.GLOBAL_RNG)
        X1 = rand(rng, n)
        X2 = rand(rng, n)
        df = DataFrame(X1 = X1, X2 = X2,
                       y = σ.(20(f.(X1) .- X2)) .> rand(rng, n))
        coerce!(df, :y => Multiclass)
    end
end

classification_data = classification_data_generator(n = 400, rng = MersenneTwister(8))

polysymbol(x, d) = Symbol(d == 0 ? "" : d == 1 ? "$x" : "$x^$d")
function poly2(data, degree)
    res = DataFrame([data.X1 .^ d1 .* data.X2 .^ d2
                     for d1 in 0:degree, d2 in 0:degree if 0 < d1 + d2 ≤ degree],
                    [Symbol(polysymbol("x₁", d1), polysymbol("x₂", d2))
                     for d1 in 0:degree, d2 in 0:degree if 0 < d1 + d2 ≤ degree])
    if hasproperty(data, :y)
        res.y = data.y
    end
    res
end;

m4 = machine(LogisticClassifier(penalty = :none),
             select(poly2(classification_data, 4), Not(:y)),
             classification_data.y) |> fit!;

m14 = machine(KNNClassifier(K = 4),
             select(classification_data, Not(:y)),
             classification_data.y) |> fit!;

path = tempname()
original = joinpath(@__DIR__, "..", "index.jl")
# so that we don't overwrite the file:
Pluto.readwrite(original, path)

@info "Loading notebook"
nb = Pluto.load_notebook(Pluto.tamepath(path));
session.notebooks[nb.notebook_id] = nb;

@info "Running notebook"
Pluto.update_save_run!(session, nb, nb.cells; run_async=false, prerender_text=true)

# nice! we ran the notebook, so we already precompiled a lot

@info "Starting HTTP server"
# next, we'll run the HTTP server which needs a bit of nasty code
t = @async Pluto.run(session)

sleep(15)
download("http://localhost:40404/")

# this is async because it blocks for some reason
@async Base.throwto(t, InterruptException())
sleep(15) # i am pulling these numbers out of thin air

@info "Warmup done"

max_stars_repo_name: ablaom/MLCourse

max_issues_repo_name: ablaom/MLCourse

max_forks_repo_name: ablaom/MLCourse


text:
struct IntervalScalingFunctionBasis
    left::Vector{Vector{Float64}}
    interior::Vector{Float64}
    right::Vector{Vector{Float64}}

    vanishing_moments::Int64
    scale::Int64
    resolution::Int64

    left_boundary::DyadicRational
    right_boundary::DyadicRational
    support::Vector{DyadicRational}

    function IntervalScalingFunctionBasis(left, interior, right, p, scale, R, left_boundary, right_boundary)
        if right_boundary <= left_boundary
            throw(error("Reconstruction interval is degenerate"))
        end

        if R < scale
            throw(error("Scale of basis functions must be greater than the R"))
        end

        length_of_support = 2*p - 1
        if length_of_support > 2^scale
            throw(DomainError(scale, "Number of vanishing moments is too large for the scale"))
        end

        left_lengths = map(length, left)
        right_lengths = map(length, right)
        if !all(length(interior) .== left_lengths .== right_lengths)
            throw(error("All basis functions must have the same length"))
        end

        if R < resolution(left_boundary) || R < resolution(right_boundary)
            throw(DomainError(R, "The R of basis functions should be at least that of the endpoints"))
        end

        s = all_dyadic_rationals(0, 1, R)
        s_translated = (right_boundary - left_boundary) * s .+ left_boundary

        new(left, interior, right, p, scale, R, left_boundary, right_boundary, s_translated)
    end
end

vanishing_moments(B::IntervalScalingFunctionBasis) = B.vanishing_moments
resolution(B::IntervalScalingFunctionBasis) = B.resolution
Base.size(B::IntervalScalingFunctionBasis) = 2^B.scale
n_eval(B::IntervalScalingFunctionBasis) = length(B.interior)

"""
    interval_scaling_function_basis(p, J, R; a, b)

Return a basis of scaling functions with `p` vanishing moments at scale `J` evaluated at resolution `R` on the interval `[a, b]`.
"""
function interval_scaling_function_basis(p::Integer, J::Integer, R::Integer; left_boundary = DyadicRational(0, 0), right_boundary = DyadicRational(1, 0))
    h = interior_filter(p)
    phi = interior_scaling_function(h, R - J)
    interior_values = 2.0^(J/2) * phi.(support(phi))

    # Ensure that the boundary values are computed
    resolution = max(R - J, ceil(Int, log2(p)) + 1)

    left = boundary_scaling_functions(LEFT, p, resolution)
    left_x = all_dyadic_rationals(support_boundaries(left[p - 1])..., R - J)
    left_values = [2.0^(J/2) * left[k].(left_x) for k in 0:p - 1]

    right = boundary_scaling_functions(RIGHT, p, resolution)
    right_x = all_dyadic_rationals(support_boundaries(right[p - 1])..., R - J)
    right_values = [2.0^(J/2) * right[k].(right_x) for k in 0:p - 1]

    IntervalScalingFunctionBasis(left_values, interior_values, right_values, p, J, R, left_boundary, right_boundary)
end

function Base.show(io::IO, B::IntervalScalingFunctionBasis)
    print(io, "Basis of scaling functions with $(vanishing_moments(B)) vanishing moments for [$(B.left_boundary), $(B.right_boundary)] at resolution $(resolution(B))")
end

function Base.getindex(B::IntervalScalingFunctionBasis, k::Integer)
    N = size(B)
    p = vanishing_moments(B)
    R = resolution(B)
    J = B.scale

    y = Vector{Float64}(undef, n_eval(B))
    x_index = evaluate!(y, B, k)

    if 1 <= k < p
        y = y[1:length(x_index)]
    elseif N - p + 1 < k <= N
        start_index = length(y) - length(x_index) + 1
        y = y[start_index:end]
    end

    return B.support[x_index], y
end

function evaluate!(y::Vector{Float64}, B::IntervalScalingFunctionBasis, k::Integer)
    N = size(B)
    p = vanishing_moments(B)
    R = resolution(B)
    J = B.scale

    if 1 <= k <= p
        x_index = 1:(2p - 1)*2^(R - J) + 1
        copy!(y, B.left[k])
    elseif p < k <= N - p
        x_index = (k - p)*2^(R - J) + 1:(p + k - 1)*2^(R - J) + 1
        copy!(y, B.interior)
    elseif N - p < k <= N
        x_index = (N - 2p + 1)*2^(R - J) + 1:2^R + 1
        copy!(y, B.right[N - k + 1])
    else
        throw(DomainError(k, "Basis functions are indexed from 1 through $N"))
    end

    return x_index
end

function reconstruct(B::IntervalScalingFunctionBasis, coef::AbstractVector)
    if length(coef) != size(B)
        throw(DimensionMismatch("The number of coefficients does not match the number of basis functions"))
    end

    x = all_dyadic_rationals(0, 1, resolution(B))
    x_translated = (B.right_boundary - B.left_boundary) * x .+ B.left_boundary

    y = Vector{Float64}(undef, n_eval(B))
    reconstruction = zeros(Float64, length(x))

    for (index, value) in enumerate(coef)
        x_index = evaluate!(y, B, index)
        r = @view reconstruction[x_index]
        LinearAlgebra.BLAS.axpy!(value, y, r)
    end

    return x_translated, reconstruction
end

function reconstruct(B::IntervalScalingFunctionBasis, coef::AbstractMatrix)
    if !all(size(coef) .== size(B))
        throw(DimensionMismatch("The number of coefficients does not match the number of basis functions"))
    end

    N = n_eval(B)
    row_values = Vector{Float64}(undef, n_eval(B))
    col_values = similar(row_values)

    x = all_dyadic_rationals(0, 1, resolution(B))
    reconstruction = zeros(Float64, length(x), length(x))

    for row_count in 1:size(coef, 1)
        row_index = evaluate!(row_values, B, row_count)
        for col_count in 1:size(coef, 2)
            col_index = evaluate!(col_values, B, col_count)

            # A loop use *a lot* less memory than a 2D BLAS.axpy!
            alpha = coef[col_count, row_count]
            r = @view reconstruction[col_index, row_index]
            for i in 1:N, j in 1:N
                r[i, j] += alpha * col_values[i] * row_values[j]
            end
        end
    end

    return reconstruction
end

@recipe function f(B::IntervalScalingFunctionBasis)
    for idx in 1:size(B)
        @series begin
            x, y = B[idx]
            float.(x), y
        end
    end
end


text:
#cd("C:/Users/LabOpto/Documents/SmartGit Projects/src")
# using Debug
# using PyPlot
using Gadfly
using DataFrames

include("../src/CDtypes.jl")
include("../src/initializations.jl")
include("../src/CDmainAlg.jl")
include("../src/GetResult.jl")
include("../src/GramMatrix.jl")
include("../src/InnerProducts.jl")
include("../src/Normalization.jl")
include("../src/utils.jl")

# get wind power data
y=readtable("wind_power.csv", header=false)
y = convert(Array , y[:,1])
y_original = copy(y)

# frequencies
f = 2*pi./collect(6:48)

@time BCD, y_path, β_best, y_best, λ_best, γ_best = l1_adaptive_trend_filter(
  y, [1, 4, 5], numλ=200, numγ=5, f=f, verbose=true
  )

# components
step = β_best[1]
# spike = β_best[2]
# slope = β_best[3]
seno = β_best[4]
cosseno = β_best[5]

# plot (Gadfly)
draw(
  SVG("step.svg", 14inch, 8inch),
  plot(
    x=1:length(step), y=step, Geom.point, Geom.line,
    Guide.xlabel("Hours"), Guide.ylabel("Step size"),
    Coord.Cartesian(xmin=0,xmax=length(y_original))
    ))
# draw(
#   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/spike.svg", 14inch, 8inch),
#   plot(x=1:length(spike),y=spike, Geom.point, Geom.line, Guide.xlabel("Hours")))
# draw(
#   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/slope.svg", 14inch, 8inch),
#   plot(x=1:length(slope),y=slope, Geom.point, Geom.line)
#   )
draw(
  SVG("trigonometric.svg", 7inch, 4.5inch),
  plot(
    layer(x=1:length(seno),y=seno, Geom.line, Theme(default_color=color("red"))),
    layer(x=1:length(cosseno), y=cosseno, Geom.line, Theme(default_color=color("blue"))),
    Guide.xlabel("frequency (ω)"), Guide.ylabel("Amplitude"), Guide.title("Trigonometric components"),
    Coord.Cartesian(xmin=0,xmax=length(f))
    ))
draw(
  SVG("fit.svg", 7inch, 4.5inch),
  plot(
    layer(x=1:length(y_original), y=y_original, Geom.line, Theme(default_color=color("red"))),
    layer(x=1:length(y_best), y=y_best, Geom.line, Theme(default_color=color("blue"))),
    Guide.xlabel("Hours"), Guide.ylabel("MW"), Guide.title("Wind Power Fit"),
    Coord.Cartesian(xmin=0,xmax=340)
    ))

# write results
writecsv("y_best.csv", y_best)
writecsv("cos.csv", cosseno)
writecsv("seno.csv", seno)
writecsv("step.csv", step)

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
export gatecount

gatecount(blk::AbstractBlock) = gatecount!(blk, Dict{Type{<:AbstractBlock}, Int}())
function gatecount!(c::Union{ChainBlock, KronBlock, PauliString, PutBlock, Bag, Concentrator, Sequence, CachedBlock}, storage::AbstractDict)
    (gatecount!.(c |> subblocks, Ref(storage)); storage)
end

function gatecount!(c::RepeatedBlock, storage::AbstractDict)
    k = typeof(content(c))
    n = length(c.locs)
    if haskey(storage, k)
        storage[k] += n
    else
        storage[k] = n
    end
    storage
end

function gatecount!(c::Union{PrimitiveBlock, Daggered, ControlBlock, ConditionBlock}, storage::AbstractDict)
    k = typeof(c)
    if haskey(storage, k)
        storage[k] += 1
    else
        storage[k] = 1
    end
    storage
end

gatecount!(c::TrivilGate, storage::AbstractDict) = storage


text:
#using colorpalettes
function colorscheme!(ec::EChart, palette::Tuple{AbstractString,Any}; reversePalette::Bool = false)

    reversePalette == true ?
		ec.color = reverse(colorpalettes[palette[1]]["$(palette[2])"]) :
			ec.color  = colorpalettes[palette[1]]["$(palette[2])"]

    return ec

end

#This is for string literals and gradients, so reversePalette doesn't make sense as argument
function colorscheme!(ec::EChart, palette::Union{AbstractString, JSFunction})

    ec.color = [palette]

    return ec

end

#For passing vector of colors
function colorscheme!(ec::EChart, palette::AbstractVector; reversePalette::Bool = false)

    reversePalette == true ? ec.color  = reverse(collect(palette)) : ec.color  = collect(palette)

    return ec

end


text:
# Fizzbuzz in Julia

for i in 1:100
    if (i % 15) == 0
        println("FizzBuzz")
    elseif (i % 3) == 0
        println("Fizz")
    elseif (i % 5) == 0
        println("Buzz")
    else
        println(i)
    end
end

max_stars_repo_name: quinnjr/fizzbuzz

max_issues_repo_name: quinnjr/fizzbuzz

max_forks_repo_name: quinnjr/fizzbuzz


text:
@testset "methods" begin
# Data
fn_phase = "data/small/Phase.nii"
fn_mag = "data/small/Mag.nii"
phase_nii = readphase(fn_phase)
mag_nii = readmag(fn_mag)
I = mag_nii .* exp.(1im .* phase_nii)
TEs = 4:4:12

# homodyne
h1 = homodyne(mag_nii, phase_nii)
h2 = homodyne(Float32.(mag_nii), Float32.(phase_nii))
h3 = homodyne(I)
h4 = homodyne(mag_nii, phase_nii; σ=[5,5])
h5 = homodyne(mag_nii, phase_nii; dims=1:3)

@test h1 == h2
@test h1 == h3
@test h1 != h4
@test h1 != h5
@test h4 != h5

# inplace test
I2 = copy(I)
@test I == I2
homodyne!(I2)
@test I != I2
@test h3 == I2

# calculateB0
B0 = calculateB0_unwrapped(romeo(phase_nii; TEs=TEs), mag_nii, TEs)
@test all(isfinite.(B0))

# make border of image noise
phase = Float32.(phase_nii)
phase[1:3,:,:,:] .= 2π .* rand.() .- π
phase[(end-2):end,:,:,:] .= 2π .* rand.() .- π
phase[:,1:3,:,:] .= 2π .* rand.() .- π
phase[:,(end-2):end,:,:] .= 2π .* rand.() .- π
phase[:,:,1:3,:] .= 2π .* rand.() .- π
phase[:,:,(end-2):end,:] .= 2π .* rand.() .- π

# getvoxelquality
vq = romeovoxelquality(phase; TEs=TEs)
@test all(isfinite.(vq))
@test size(vq) == size(phase_nii)[1:3]

# mask_from_voxelquality
mask = mask_from_voxelquality(vq, 0.5)
@test mask isa AbstractArray{<:Bool}
@test !all(mask .== true)
@test !all(mask .== false)

end


text:
#module HierarchicalBetaBinomial

import Distributions: Beta, Binomial, Normal
import LinearAlgebra

include("./src/types.jl")
include("./src/MH.jl")

#end


text:
# Methods where dims are an argument. This is an experimental approach
# that might need to change.
# targeting underscore _methods so we can use dispatch on the dims arg

# TODO hanbdle rebuild in the array dispatch, the dimensions are probably wrong in some cases.
for (mod, fname) in ((:Base, :sum), (:Base, :prod), (:Base, :maximum), (:Base, :minimum), (:Statistics, :mean))
    _fname = Symbol('_', fname)
    @eval begin
        @inline ($mod.$_fname)(a::AbstractArray{T,N}, dims::AllDimensions) where {T,N} =
            ($mod.$_fname)(a, dimnum(a, dims))
        @inline ($mod.$_fname)(f, a::AbstractArray{T,N}, dims::AllDimensions) where {T,N} =
            ($mod.$_fname)(f, a, dimnum(a, dims))
    end
end

for fname in (:std, :var)
    _fname = Symbol('_', fname)
    @eval begin
        @inline (Statistics.$_fname)(a::AbstractArray, corrected::Bool, mean, dims::AllDimensions) =
            (Statistics.$_fname)(a, corrected, mean, dimnum(a, dims))
    end
end

Statistics._median(a::AbstractArray, dims::AllDimensions) =
    Statistics._median(a, dimnum(a, dims))
Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractArray, dims::AllDimensions) =
    Base._mapreduce_dim(f, op, nt, A, dimnum(A, dims))
# Unfortunately Base/accumulate.jl kwargs methods all force dims to be Integer.
# accumulate wont work unless that is relaxed, or we copy half of the file here.
Base._accumulate!(op, B, A, dims::AllDimensions, init::Union{Nothing, Some}) =
    Base._accumulate!(op, B, A, dimnum(A, dims), init)

Base._dropdims(a::AbstractArray, dim::Union{AbDim,Type{<:AbDim}}) =
    rebuildsliced(a, Base._dropdims(a, dimnum(a, dim)), dims2indices(a, basetype(dim)(1)))
Base._dropdims(a::AbstractArray, dims::AbDimTuple) =
    rebuildsliced(a, Base._dropdims(a, dimnum(a, dims)),
                  dims2indices(a, Tuple((basetype(d)(1) for d in dims))))

# TODO cov, cor mapslices, eachslice, reverse, sort and sort! need _methods without kwargs in base so
# we can dispatch on dims. Instead we dispatch on array type for now, which means
# these aren't usefull unless you inherit from AbDimArray.

Base.mapslices(f, a::AbDimArray; dims=1, kwargs...) = begin
    dimnums = dimnum(a, dims)
    data = mapslices(f, parent(a); dims=dimnums, kwargs...)
    rebuildsliced(a, data, dims2indices(a, reducedims(DimensionalData.dims(a, dimnums))))
end

# This is copied from base as we can't efficiently wrap this function
# through the kwarg with a rebuild in the generator. Doing it this way
# wierdly makes it 2x faster to use a dim than an integer.
if VERSION > v"1.1-"
    Base.eachslice(A::AbDimArray; dims=1, kwargs...) = begin
        if dims isa Tuple && length(dims) == 1
            throw(ArgumentError("only single dimensions are supported"))
        end
        dim = first(dimnum(A, dims))
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
        return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
    end
end

for fname in (:cor, :cov)
    @eval Statistics.$fname(a::AbDimArray{T,2}; dims=1, kwargs...) where T = begin
        newdata = Statistics.$fname(parent(a); dims=dimnum(a, dims), kwargs...)
        I = dims2indices(a, dims, 1)
        newdims, newrefdims = slicedims(a, I)
        rebuild(a, newdata, (newdims[1], newdims[1]), newrefdims)
    end
end

@inline Base.reverse(a::AbDimArray{T,N}; dims=1) where {T,N} = begin
    dnum = dimnum(a, dims)
    # Reverse the dimension. TODO: make this type stable
    newdims = revdims(DimensionalData.dims(a), dnum)
    # Reverse the data
    newdata = reverse(parent(a); dims=dnum)
    rebuild(a, newdata, newdims, refdims(a))
end

@inline revdims(dimstorev::Tuple, dnum) = begin
    dim = dimstorev[end]
    if length(dimstorev) == dnum
        dim = rebuild(dim, reverse(val(dim)))
    end
    (revdims(Base.front(dimstorev), dnum)..., dim)
end
@inline revdims(dims::Tuple{}, i) = ()

for fname in [:permutedims, :transpose, :adjoint]
    @eval begin
        @inline Base.$fname(a::AbDimArray{T,2}) where T =
            rebuild(a, $fname(parent(a)), reverse(dims(a)), refdims(a))
    end
end

Base.permutedims(a::AbDimArray{T,N}, perm) where {T,N} =
    rebuild(a, permutedims(parent(a), dimnum(a, perm)), permutedims(dims(a), perm))


text:
using JuMP
using Ipopt
using LinearAlgebra

X = randn(10,1000)
Xn = sqrt.(sum(X.^2, dims=2))
X ./= Xn
ϵ = 1e-1
λ = 1.0
α = 1.0
n_eigs = 5
v, V = eigen(X'*X + ϵ*I)
Γ = V[:, (end-n_eigs+1):end]

μ = rand(10)
U = copy(Γ)

function Ml(μ, X, ϵ)
    M = ϵ*I + X'*Diagonal(μ)*X
    Mc = cholesky(Hermitian(M))
    return (M, Mc)
end

Ml(μ) = Ml(μ, X, ϵ)

_, Mc = Ml(μ)

f(U) = α*sum((U.-Γ).^2) + λ*tr(U'*(Mc\U))

function g!(gvec, U, m, n)
    for i = 1:n
        gvec[((i-1)*m+1):i*m] .= 2*α*(U[1:m, i].-Γ[1:m, i]) .+ 2*λ*(Mc\U[1:m, i])
    end
end

m, n = size(U)

L = LinearIndices((m, n))
jump_wrapper_f(x...) = f([x[L[i,j]] for i in 1:m, j in 1:n])
jump_wrapper_g!(gvec, x...) = g!(gvec, [x[L[i,j]] for i in 1:m, j in 1:n], m, n)

sumprod(x...) = sum(x[i]*x[i+m] for i in 1:m)
function sumprod_grad!(gvec, x...)
    for i in 1:m
        gvec[i] = x[i+m]
        gvec[i+m] = x[i]
    end
end

model = Model(optimizer_with_attributes(Ipopt.Optimizer))
JuMP.register(model, :f, m*n, jump_wrapper_f, jump_wrapper_g!)
JuMP.register(model, :sumprod, 2*m, sumprod, sumprod_grad!)

@variable(model, x[1:m, 1:n])

set_start_value.(x, U)

# This is the problem section because of the indexing expressions being splatted

for i = 1:n
    for j = i:n
        x_i = x[1:m, i]
        x_j = x[1:m, j]
        @NLconstraint(model, sumprod(x_i..., x_j...) == ((i==j) ? 1 : 0))
    end
end

@NLobjective(model, Min, f(x...))

optimize!(model)

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
# To test a simplified quadrotor.
include("../src/plotSol.jl")
include("../src/simulate.jl")
# include("../src/linearize.jl")
include("../src/nRotor.jl")
# include("../src/trim_kronLazy.jl");
# include("../src/trim_kron.jl")
# include("../src/trim_FinDiff.jl")

# using nRotor.jl, generate quadrotor joint tree
qRotor_j = gen_nRotor(4);
g =([0.0,0.0,-9.81]); # Gravity Vector.

## Simulate
tEnd = 0.1;
tSpan = 0.01;
tSim, solFinal = simulate(tEnd, tSpan, qRotor_j, g=g);

solQuad = solFinal[1];
solProp1 = solFinal[2];
solProp2 = solFinal[3];
solProp3 = solFinal[4];
solProp4 = solFinal[5];

##
ωCube, ωProp1, ωRel1, jointLoc1 = checkRevJoint(solQuad, solProp1, qRotor_j[2].pos1, qRotor_j[2].pos2);
ωCube, ωProp2, ωRel2, jointLoc2 = checkRevJoint(solQuad, solProp2, qRotor_j[3].pos1, qRotor_j[3].pos2);
ωCube, ωProp3, ωRel3, jointLoc3 = checkRevJoint(solQuad, solProp3, qRotor_j[4].pos1, qRotor_j[4].pos2);
ωCube, ωProp4, ωRel4, jointLoc4 = checkRevJoint(solQuad, solProp4, qRotor_j[5].pos1, qRotor_j[5].pos2);

## plots
# errors
# close("all");
# plotErrNorm(tSim,solQuad.β)
# plotErrNorm(tSim,solProp1.β)
# plotErrNorm(tSim,solProp2.β)
# plotErrNorm(tSim,solProp3.β)
# plotErrNorm(tSim,solProp4.β)
# plotPos(tSim,jointLoc1);
# plotPos(tSim,jointLoc2);
# plotPos(tSim,jointLoc3);
# plotPos(tSim,jointLoc4);
#
# # trajectory
# plotPos(tSim, solQuad.r, "Quad");
# plotPos(tSim, solProp1.r, "Prop1");
# plotPos(tSim, solProp2.r, "Prop2");
# plotPos(tSim, solProp3.r, "Prop3");
# plotPos(tSim, solProp4.r, "Prop4");
#
# # angular velocities
# plotAngVel(tSim, ωCube)
# plotAngVel(tSim, ωProp1)

## find trim point and obtain linear plant
x0, u0 = getXU_0(qRotor_j);
nB = qRotor_j[end].RB2.bodyID;
ix= zeros(Integer,20*(nB-1)+1);
# WITH no fault in actuators, hover condition
ix[14*(nB-1)+1:20*(nB-1)] .= 1;
trim_U = u0;
iy = zeros(Integer,7*(nB-1));

out = trim_kronLazy(qRotor_j,g,ix=ix, iy=iy);
# out = trim_FinDiff(qRotor_j,g,ix=ix, iy=iy);
trim_x, trim_u, gam = out;
println("trim_x = ", trim_x)
trim_X = [qRotor_j[1].RB1.x;trim_x]
# trim_U = [zeros(6) reshape(trim_u,(6,nB-1))];
# println("trim_U = ", trim_U[:,2])
# trim_U = zeros(6,nB);
println("qconstr = ", norm(trim_x[4:7])-1)
println("ẍ =", norm(fxdot(trim_X,trim_U, qRotor_j, g)))

## old
## find linear system
x0, u0 = getXU_0(qRotor_j);
x0_woIn = x0[15:end]; u0_woIn = u0[:,2:end]; # without the inertial frame
println()
println("####Linearizing####")
A, Bu = linearize(x0, u0, qRotor_j, g);
Cy = Matrix{Float64}(I,size(A)...); # full state measurement
Cz = Matrix{Float64}(I,size(A)...); # full state output
Dzu = zeros(size(Bu));
## lti system (nominal, no disturbance or faults)
using PyCall
pyc = pyimport("control");
Usys = zeros(length(u0_woIn),length(tSim));
for i=1:length(tSim)
  Usys[:,i] = u0_woIn[:];
end
sysNom = pyc.ss(A, Bu, Cz, Dzu); # nominal system at hover
T, yout, xout = pyc.forced_response(sysNom, permutedims(tSim), Usys, x0_woIn);

## with disturbance
using Random
rng = MersenneTwister(1234);
ξ = randn(1,length(tSim));
Uξsys = [Usys; ξ];

Dyξ = ones(length(x0_woIn)); # disturbance in all states
# no input in measurement
Bξ = ones(length(x0_woIn)); # disturbance in all states
Buξ = [Bu Bξ];
Dzξ = ones(length(x0_woIn)); # disturbance in all states
Dzuξ = [Dzu Dzξ];
sysDist = pyc.ss(A, Buξ, Cz, Dzuξ); # nominal system at hover
T2, yout2, xout2 = pyc.forced_response(sysDist, permutedims(tSim), Uξsys, x0_woIn);

## with faults and disturbance
flt_act = ones(length(u0_woIn)); # actuator fault vector; 1 indicates no fault
act_pos = 9:6:27; # position of thrust in the u vector
flt_act[act_pos] = [1.0;0.5;1.0;0.5]; # partial fault in rotors 2 and 4
Σa = diagm(flt_act);
# BuΣa = Bu*Σa;
Bfξ = [Bu*Σa Bξ];
Dfξ = [Dzu*Σa Dzξ];

flt_sen = ones(size(Cy,1)); # sensor fault vector
Σs = diagm(flt_sen);
Σs_Cy = Σs*Cy; Σs_Dyξ = Σs*Dyξ;

sysFlt = pyc.ss(A, Bfξ, Cz, Dfξ); # system with fault
Tf, yout_f, xout_f = pyc.forced_response(sysFlt, permutedims(tSim), Uξsys, x0_woIn);

## old
# clearconsole()
# j1 = Nothing
# j2 = Nothing
# j3 = Nothing
# j4 = Nothing
# j5 = Nothing
# ##
# mb = 0.155; mp = 0.02;
#
# Ib =  [ 0.0009    0.0000   -0.0000;
#         0.0000    0.0053    0.0000;
#        -0.0000    0.0000    0.0054]#*(10^-3)
#
# Ip =    [ 0.0029    0.0000    0.0000;
#           0.0000    0.5824    0.0000;
#           0.0000    0.0000    0.5805]*(10^-4)
#
# InFrame = InertialFrameAsRB()
# QuadCube = RigidBody(mb,Ib,2)
# Prop1 = RigidBody(mp,Ip,3)
# Prop2 = RigidBody(mp,Ip,4)
# Prop3 = RigidBody(mp,Ip,5)
# Prop4 = RigidBody(mp,Ip,6)
#
# x0Cube = [zeros(3);[1;zeros(3)];zeros(3);zeros(4)]
# initialiseRigidBody!(QuadCube,x0Cube)
#
# x0Prop1 = [[0.5;-0.5;0.1];[1;zeros(3)];zeros(3);zeros(4)]
# initialiseRigidBody!(Prop1,x0Prop1)
# x0Prop2 = [[0.5;0.5;0.1];[1;zeros(3)];zeros(3);zeros(4)]
# initialiseRigidBody!(Prop2,x0Prop2)
# x0Prop3 = [[-0.5;-0.5;0.1];[1;zeros(3)];zeros(3);zeros(4)]
# initialiseRigidBody!(Prop3,x0Prop3)
# x0Prop4 = [[-0.5;0.5;0.1];[1;zeros(3)];zeros(3);zeros(4)]
# initialiseRigidBody!(Prop4,x0Prop4)
#
# axis = [0.0 0.0 1.0][:]
# rjCube1 = [0.5 -0.5 0.1][:]
# rjCube2 = [0.5 0.5 0.1][:]
# rjCube3 = [-0.5 -0.5 0.1][:]
# rjCube4 = [-0.5 0.5 0.1][:]
# rjProp = [0.0 0.0 0.0][:];
#
# j1 = Joint(InFrame,QuadCube,zeros(3),zeros(3));
# j2 = Joint(QuadCube,Prop1,rjCube1,rjProp,type = "Revolute",axis = axis);#,jointTorque = [0.0;0.0;0.0001]);
# j3 = Joint(QuadCube,Prop2,rjCube2,rjProp,type = "Revolute",axis = axis);
# j4 = Joint(QuadCube,Prop3,rjCube3,rjProp,type = "Revolute",axis = axis);
# j5 = Joint(QuadCube,Prop4,rjCube4,rjProp,type = "Revolute",axis = axis);
#
# j = ([j1,j2,j3,j4,j5]);
# ## Linearization
# # x0, u0 = getXU_0(j);
# # A,B = linearizeDiff(x0, u0, j);
# # println("A = ",A)
# # println("B = ",B)
# # sleep(1200)
#
# ## Simulation
# tEnd = 0.1;
# tSpan = 0.01;
# g =([0.0,0.0,0.0]) # Gravity Vector.
# tSim, solFinal = simulateDiff(tEnd,tSpan,j,g=g);
#
# solQuad = solFinal[1];
# solProp1 = solFinal[2];
# solProp2 = solFinal[3];
# solProp3 = solFinal[4];
# solProp4 = solFinal[5];
#
# ##
# ωCube, ωProp1, ωRel1, jointLoc1 = checkRevJoint(solQuad, solProp1, rjCube1, rjProp);
# ωCube, ωProp2, ωRel2, jointLoc2 = checkRevJoint(solQuad, solProp2, rjCube2, rjProp);
# ωCube, ωProp3, ωRel3, jointLoc3 = checkRevJoint(solQuad, solProp3, rjCube3, rjProp);
# ωCube, ωProp4, ωRel4, jointLoc4 = checkRevJoint(solQuad, solProp4, rjCube4, rjProp);
# #
#
# # Check if joint location has moved
# # jointLoc1 = Matrix{Float64}(undef,length(tSim),3);
#
# # ωCube = Matrix{Float64}(undef,length(tSim),3);
# # ωProp1 = Matrix{Float64}(undef,length(tSim),3);
# # ωProp1InCube = Matrix{Float64}(undef,length(tSim),3);
# # ωRel = Matrix{Float64}(undef,length(tSim),3);
# # for i=1:length(tSim)
# #     QuadCube.dcm = quat2dcm(solQuad.β[i,:])
# #     ωCube[i,:] = angVel(solQuad.β[i,:],solQuad.βdot[i,:])
# #     ωProp1[i,:] = angVel(solProp1.β[i,:],solProp1.βdot[i,:])
# #     ωProp1InCube[i,:] = quat2dcm(solQuad.β[i,:])*transpose(quat2dcm(solProp1.β[i,:]))*ωProp1[i,:]
# #     ωRel[i,:] = ωProp1InCube[i,:] - ωCube[i,:]
# #     jointLoc1[i,:] = solQuad.r[i,:] + transpose(QuadCube.dcm)*rjCube1 - solProp1.r[i,:]
# # end
# ##
# close("all");
# # plotErrNorm(tSim,solQuad.β)
# # plotErrNorm(tSim,solProp1.β)
# plotPos(tSim,jointLoc2);
# # plotVel(tSim,solQuad.v  - solProp.v);
# # plotAngVel(tSim,ωCube);
# # plotAngVel(tSim,ωRel1);
# # plotAngVel(tSim,ωProp1)
# # plotQuatDot(tSim, solQuad.βdot)
# ##
# # close("all");
# # plotPos(tSim, solQuad.r, "Quad");
# # # plotQuat(tSim, solQuad.β)
# # # plotQuatDot(tSim, solQuad.βdot)
# # # plotErrNorm(tSim, solQuad.β)
# # plotPos(tSim, solProp1.r, "Prop1");
# # # plotVel(tSim, solProp1.v, "Prop1");
# # plotPos(tSim, solProp2.r, "Prop2");
# # plotPos(tSim, solProp3.r, "Prop3");
# # plotPos(tSim, solProp4.r, "Prop4");
# ##
# plotAngVel(tSim, ωCube)
# plotAngVel(tSim, ωProp1)
# plotAngVel(tSim, ωProp3)
#

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
using JSON, HTTP

const HEADERS = Dict("Authorization" => "Bearer " * ENV["FLUXBOT_BUILDKITE_SECRET"])
const ORG_SLUG = "julialang"
const BUILDKITE_BASE_URL = "https://api.buildkite.com/v2/organizations/$ORG_SLUG/pipelines/"

### Jobs API

function create_job(project_slug, commit, branch, env = Dict(); headers = HEADERS)

  data = Dict("commit" => commit, "branch" => branch, "env" => env)
  url = joinpath(BUILDKITE_BASE_URL, project_slug, "builds")
  req = HTTP.post(url, headers, JSON.json(data))
  JSON.parse(String(req.body))
end

function get_job_logs(project_slug, build_number, job_id; headers = HEADERS)
  url = joinpath(BUILDKITE_BASE_URL, project_slug, "builds", build_number, "jobs", job_id, "log")
  req = HTTP.get(url, headers)
  JSON.parse(String(req.body))
end

### Builds API

function get_build(project_slug, build_number; headers = HEADERS)
  url = joinpath(BUILDKITE_BASE_URL, project_slug, "builds", build_number)
  req = HTTP.get(url, headers = headers)
  JSON.parse(String(req.body))
end
get_build(project_slug; kwargs...) = get_build(project_slug, ""; kwargs...)
get_build(project_slug, build::Int; kw...) = get_build(project_slug, string(build); kw...)

function create_build(project_slug, commit, branch; headers = HEADERS)
  url = joinpath(BUILDKITE_BASE_URL, project_slug, "builds")
  data = Dict("commit" => commit, "branch" => branch)
  req = HTTP.post(url, headers, JSON.json(data))
  JSON.parse(String(req.body))
end

function cancel_build(project_slug, build_number; headers = HEADERS)
  url = joinpath(BUILDKITE_BASE_URL, project_slug, "builds", build_number, "cancel")
  req = HTTP.put(url, headers)
  JSON.parse(String(req.body))
end

function rebuild_build(project_slug, build_number; headers = HEADERS)
  url = joinpath(BUILDKITE_BASE_URL, project_slug, "builds", build_number, "rebuild")
  req = HTTP.put(url, headers)
  JSON.parse(String(req.body))
end