text:
"""
Normalize time-domain seismic data.

# Arguments

* `data::TD` : input data
* `attrib::Symbol` : decide kind of normalization
  * `=:recrms` the record at every receiver is normalized with its RMS value
  * `=:recmax` the record at every receiver is normalized with its maximum value

# Return

* normalized data as `TD`
"""
function TD_normalize(data::TD, attrib::Symbol=:recrms)
	datan = deepcopy(data);
	TD_normalize!(datan, attrib)
	return datan
end
function TD_normalize!(data::TD, attrib::Symbol=:recrms)
	nr = data.ageom.nr;	nss = data.ageom.nss;	nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		scs=vecnorm(dd,2)

		for ir = 1:nr[iss]
			ddv=view(dd, :, ir)

			if(attrib == :recrms)
				sc=vecnorm(ddv,2)
				rmul!(ddv,inv(sc))
			elseif(attrib == :recmax)
				sc=vecnorm(ddv,Inf)
				rmul!(ddv,inv(sc))
			elseif(attrib == :srcrms)
				rmul!(ddv,inv(scs))
			else
				error("invalid attrib")
			end
		end
	end
	return data
end

function TD_filter!(data::TD; fmin=nothing, fmax=nothing)

	if((fmin===nothing) && (fmax===nothing))
		return nothing
	end
	fs = inv(step(data.tgrid));
	designmethod = Butterworth(4);
	filtsource = Bandpass(fmin, fmax; fs=fs);

	nr = data.ageom.nr;	nss = data.ageom.nss;	nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		scs=vecnorm(dd,2)

		for ir = 1:nr[iss]
			ddv=view(dd, :, ir)

			filt!(ddv, digitalfilter(filtsource, designmethod), ddv)
		end
	end
	return data
end

"""
Construct TD using data at all the unique receiver positions
for all supersources.

* `d::Array{Float64}` : the data matrix ordered in order such that time-domain modelling schemes are fast, i.e., [irec,ifield,it,nss]

"""
function TD_urpos(d::Array{Float64},
		   fields::Vector{Symbol},
		   tgrid::StepRangeLen,
		   acq::AGeom,
		   nur::Int64,
		   urpos::Tuple{Array{Float64,1},Array{Float64,1}
		  }
		   )
	dout = [zeros(length(tgrid),acq.nr[iss]) for iss=1:acq.nss, ifield=1:length(fields)]

	for ifield=1:length(fields), iss=1:acq.nss, ir=1:acq.nr[iss]
		# find index in urpos
		irr=find([[urpos[1][i]-acq.rz[iss][ir],
		       urpos[2][i]-acq.rx[iss][ir]] == [0., 0.,] for i in 1:nur])

		dout[iss, ifield][:,ir] = d[irr[1],ifield, :,iss]
	end

	return TD(dout, fields, tgrid, acq)

end

"""
Apply source and receiver coupling functions to TD.
Currently, only source filters are applied.

# Arguments

* `s::TD` : input data
* `r::TD` : input data
* `w::Coupling.TD` : input source and receiver filters
* `attrib::Symbol` : attribute to
  * `=:s` to apply `w` to `r` and modify `s`
  * `=:r` to apply adjoint of `w` to `s` and modify `r`
  * `=:w` modify `w` using `r` and `s`

TODO: need to work on parallelization and speed up here
"""
function TDcoup!(
               s::TD,
	       r::TD,
	       w::Coupling.TD,
	       attrib::Symbol
	       )
	nr = r.ageom.nr;	nss = r.ageom.nss;	nt = length(r.tgrid);
	fields = (w.fields == r.fields == s.fields) ? w.fields : error("different fields")
	sv=zeros(length(s.tgrid))
	rv=zeros(length(r.tgrid))
	wv=zeros(length(w.tgridssf))
	for ifield = 1:length(fields), iss = 1:nss, ir = 1:nr[iss]
		# receiver coupling
	#	conv!(s.d[iss, ifield][:, ir],r.d[iss, ifield][:, ir],
#		 w.rf[iss, ifield][:,ir],:s)
		# source coupling
		sv=s.d[iss, ifield][:, ir]
		rv=r.d[iss, ifield][:, ir]
		wv=w.ssf[iss, ifield]
#		conv!(sv, rv, wv, attrib)
		s.d[iss, ifield][:, ir]=copy(sv)
		r.d[iss, ifield][:, ir]=copy(rv)
		w.ssf[iss, ifield][:]=copy(wv)
	end
end # TDcoup


text:
using StaticArrays
using IntervalArithmetic
import IntervalArithmetic: interval

using LinearAlgebra

import Base: +, *, ^, -, sqrt, inv

"""
Affine form with center `c`, affine components `γ` and error `Δ`.

Variant where Δ is an interval
"""
struct Affine{N,T<:AbstractFloat}
    affine::Aff{N,T}
    range::Interval{T}
end

interval(x::Affine) = x.range
range(x::Affine) = x.range

eltype(::Affine{N, T}) where {N, T} = T
zero(::Affine{N, T}) where {N, T} = Affine(Interval(zero(T)))
zero(::Type{Affine{T}}) where {N, T} = Affine(Interval(zero(T)))

one(::Affine{T}) where T = Affine(Interval(one(T)))
one(::Type{Affine{T}}) where T = Affine(Interval(one(T)))

function Base.show(io::IO, C::Affine{N,T}) where {N,T}
    print(io, "affine=", C.affine, "; range=", C.range)
end

==(C::Affine, D::Affine) = C.affine == D.affine && C.range == D.range

"""
Make an `Affine` based on an interval, which is number `i` of `n` total variables.
"""
function Affine(X::Interval{T}, n, i) where {T}
    return Affine(Aff(X, n, i), X)
end

Affine(X::Interval) = Affine(X, 1, 1)
Affine(X::Number) = Affine(Interval(X), 1, 1)

affine(Xs::Interval...) = Affine.(Xs, length(Xs), 1:length(Xs))

for op in (:+, :*, :-)
    @eval function $op(x::Affine, y::Affine)
        affine = $op(x.affine, y.affine)

        range = $op(x.range, y.range)

        range = range ∩ interval(affine)

        return Affine(affine, range)
    end
end

for op in (:sqrt, :inv)
    @eval function $op(x::Affine)
        affine = $op(x.affine, x.range)

        range = $op(x.range)

        range = range ∩ interval(affine)

        return Affine(affine, range)
    end
end

*(x::Affine, α::Real) = Affine(α*x.affine, α*x.range)
*(α::Real, x::Affine) = x * α

+(x::Affine, α::Real) = Affine(α+x.affine, α + x.range)
+(α::Real, x::Affine) = x + α

-(x::Affine) = Affine(-x.affine, -x.range)
-(x::Affine, α::Real) = Affine(x.affine - α, x.range - α)
-(α::Real, x::Affine) = α + (-x)

/(x::Affine, α::Real) = Affine(x.affine/α, x.range/α)

function ^(x::Affine, n::Integer)

    invert = false

    if n < 0
        invert = true
        n = -n
        # @show n
    end

    result = Base.power_by_squaring(x, n)

    if invert
        result = inv(result)
    end

    return result
end
Base.literal_pow(::typeof(^), x::Affine, ::Val{p}) where {T,p} = x^p

ext: jl

max_stars_repo_name: dpsanders/AffineArithmetic

max_issues_repo_name: dpsanders/AffineArithmetic

max_forks_repo_name: dpsanders/AffineArithmetic


text:
@info "Análise Estrutural"
# Matriz de rigidez dos elementos
k1 = [300  -300
    -300   300]

k2 = [
    600  -600
    -600   600
]

k3 = [
    800  -800
    -800   800
]

k4 = [
    500  -500
    -500   500
]

k5 = [
    600  -600
    -600   600
]

k6 = [
    800  -800
    -800   800
]

# Vetor de forças aplicadas nos Nós
f =[0, 500, 500, 600, 1000, 0]
# Vetor de deslocamentos nodais
u = [0, 0, 0, 0, 0, 2]

# Nós
A = 1
B = 2
C = 3
D = 4
E = 5
F = 6

# Relacionamento entre elementos e nós
ak1 = [[A,A], [A,B], [B,A], [B,B]]
ak2 = [[B,B], [B,C], [C,B], [C,C]]
ak3 = [[B,B], [B,D], [D,B], [D,D]]
ak4 = [[B,B], [B,E], [E,B], [E,E]]
ak5 = [[D,D], [F,D], [D,F], [F,F]]
ak6 = [[E,E], [E,F], [F,E], [F,F]]

# Matriz de rigidez global
k = zeros(6, 6)

# Mapeia a matriz de rigidez local com a matriz de rigidez global
function mapear(a, b)
    l = 1
    for i in 1:2
        for j in 1:2
            k[a[l][1,1],a[l][2,1]] += b[i, j]
            l+=1
        end
    end
end

# Monta matriz de rigidez global
mapear(ak1, k1)
mapear(ak2, k2)
mapear(ak3, k3)
mapear(ak4, k4)
mapear(ak5, k5)
mapear(ak6, k6)

# Resolve sistema de equações para achar os deslocamentos nodais
# é retirada da matriz a linha e coluna correspondente ao nó que não tem deslocamento
@show _k = (k[1:end .!=A, 1:end .!=A])[1:end .!=E, 1:end .!=E] \ (f[1:end .!=A])[1:end .!=E]

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
using JuliaTutor
using Test

@testset "JuliaTutor.jl" begin
    # There is almost nothing to actually test in this package...
    # It will be interesting to see what people come up with.
end


text:
module Mos_Model
using JuMP
using KNITRO

function mos_model(glc_ext)
    model = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "honorbnds" => 1,
            "ms_maxsolves" => 10,
        ), # might need more resolves when model gets more complicated
    )

    # make functions
    ϵ = 1e-6
    fmax(x) = (x + sqrt(ϵ + x^2)) / 2 # max(0, x) approximation (this creates a normal Julia function)
    register(model, :fmax, 1, fmax; autodiff = true) # JuMP emits a warning telling you to do this
    fmin(x) = -fmax(-x) # min(0, x) approximation
    register(model, :fmin, 1, fmin; autodiff = true) # JuMP emits a warning telling you to do this

    # bounds
    lb_enz = 0.0 # [g enz / g DW cell]
    ub_enz = 1.0 # [g enz / g DW cell]
    lb_met = log(1e-6) # [log(1 uM)]
    ub_met = log(10e-3) # [log(10 mM)]
    lb_dg = -100 # [kJ/mol]
    ub_dg = 100 # [kJ/mol]
    lb_v = -50 # [mmol/gDW/h]
    ub_v = 50 # [mmol/gDW/h]

    @variables model begin
        # biomass [1/h] (g6p + ATP -> ADP)
        0 <= mu <= 10

        # Enzyme concentrations [g enz / g DW]
        lb_enz <= pts <= ub_enz # enzyme for: glc external + ATP -> ADP + g6p
        lb_enz <= emp <= ub_enz # enzyme for: g6p + 2NAD + 2ADP -> 2NADH + 2ATP + 2pep
        lb_enz <= pyk <= ub_enz # enzyme for: pep + ADP -> ATP + pyr
        lb_enz <= ldh <= ub_enz # enzyme for: pyr + NADH -> NAD + lac
        lb_enz <= burn <= ub_enz # enzyme for: ATP -> ADP
        lb_enz <= lp <= ub_enz # enzyme for: lac internal -> lac external (lactate permease)

        # Concentrations of internal metabolites [log(M)]
        lb_met <= g6p <= ub_met # glucose-6-phosphate
        lb_met <= pyr <= ub_met # pyruvate
        lb_met <= pep <= ub_met # phosphoenolpyrovate
        lb_met <= lac <= ub_met # lactate
        lb_met <= atp <= ub_met # ATP
        lb_met <= adp <= ub_met # ADP
        lb_met <= nad <= ub_met # NAD
        lb_met <= nadh <= ub_met # NADH

        # Fluxes [mmol/gDW/h]
        lb_v <= v_pts <= ub_v
        lb_v <= v_emp <= ub_v
        lb_v <= v_pyk <= ub_v
        lb_v <= v_ldh <= ub_v
        lb_v <= v_burn <= ub_v
        lb_v <= v_lp <= ub_v

        # Thermodynamic variables
        lb_dg <= dg_pts <= ub_dg
        lb_dg <= dg_emp <= ub_dg
        lb_dg <= dg_pyk <= ub_dg
        lb_dg <= dg_ldh <= ub_dg
        lb_dg <= dg_burn <= ub_dg
        lb_dg <= dg_lp <= ub_dg
    end

    @NLparameters model begin
        # media conditions
        glc_ext == log(50e-3) # [log(50 mM)]
        lac_ext == log(10e-6) # [log(10 uM)]

        # capacity constraint
        total_proteome_mass_fraction == 0.26 # [g enz/gDW]

        # maintenance requirement
        min_burn_flux == 1 # [mmol/gDW/h]

        # thermodynamic properties
        RT == 8.3145e-3 * 298.15 # gas constant [kJ/K/mol] * temperature [K]

        # enzyme specific activities (related to kcats) [mmol/g enz/h]
        kcat_pts == 213.75
        kcat_emp == 126
        kcat_pyk == 126
        kcat_ldh == 31
        kcat_burn == 10
        kcat_lp == 100

        # Standard Gibbs energy of reaction [kJ/mol]
        dg0_pts == -16.7
        dg0_emp == -4.71
        dg0_pyk == -31.7
        dg0_ldh == -23.7
        dg0_burn == -57
        dg0_lp == -10.3
    end

    @NLconstraints model begin
        # ΔG at concentrations
        dg_pts == dg0_pts + RT * (adp + g6p - atp - glc_ext)
        dg_emp == dg0_emp + RT * (2 * nadh + 2 * atp + 2 * pep - 2 * adp - 2 * nad - g6p)
        dg_pyk == dg0_pyk + RT * (pyr + atp - adp - pep)
        dg_ldh == dg0_ldh + RT * (lac + nad - nadh - pyr)
        dg_burn == dg0_burn + RT * (adp - atp)
        dg_lp == dg0_lp + RT * (lac_ext - lac)

        # flux bounds due to kinetics and thermodynamics
        pts * kcat_pts * fmin(exp(-dg_pts / RT) - 1) <= v_pts
        emp * kcat_emp * fmin(exp(-dg_emp / RT) - 1) <= v_emp
        pyk * kcat_pyk * fmin(exp(-dg_pyk / RT) - 1) <= v_pyk
        ldh * kcat_ldh * fmin(exp(-dg_ldh / RT) - 1) <= v_ldh
        burn * kcat_burn * fmin(exp(-dg_burn / RT) - 1) <= v_burn
        lp * kcat_lp * fmin(exp(-dg_lp / RT) - 1) <= v_lp

        v_pts <= pts * kcat_pts * fmax(1 - exp(dg_pts / RT))
        v_emp <= emp * kcat_emp * fmax(1 - exp(dg_emp / RT))
        v_pyk <= pyk * kcat_pyk * fmax(1 - exp(dg_pyk / RT))
        v_ldh <= ldh * kcat_ldh * fmax(1 - exp(dg_ldh / RT))
        v_burn <= burn * kcat_burn * fmax(1 - exp(dg_burn / RT))
        v_lp <= lp * kcat_lp * fmax(1 - exp(dg_lp / RT))

        # mass balance constraints
        v_pts - v_emp == mu # g6p
        2 * v_emp - v_pyk == 0 # pep
        v_pyk - v_ldh == 0 # pyr
        v_ldh - v_lp == 0 # lac
        2 * v_emp + v_pyk - v_pts - v_burn == mu # atp
        v_pts + v_burn - 2 * v_emp - v_pyk == -mu # adp
        v_ldh - 2 * v_emp == 0 # nad
        2 * v_emp - v_ldh == 0 # nadh

        # density constraint(s) (can add more, e.g. membrane)
        pts + emp + pyk + ldh + lp + burn <= total_proteome_mass_fraction

        # minimum maintenance
        min_burn_flux <= v_burn

        # measured ratio constraints (not strictly necessary)
        atp == log(10) + adp # atp/adp = 10, but remember that the concentration variables are logged
        nadh == log(0.2) + nad
    end

    @objective(model, Max, mu)
    optimize!(model)
    objective_value(model)

    results_mo = Dict(
        "mu" => value(mu),
        "glc_ext" => value(glc_ext),
        "g6p" => value(g6p),
        "pyr" => value(pyr),
        "pep" => value(pep),
        "lac" => value(lac),
        "atp" => value(atp),
        "adp" => value(adp),
        "nad" => value(nad),
        "nadh" => value(nadh),
        "v_pts" => value(v_pts),
        "v_emp" => value(v_emp),
        "v_pyk" => value(v_pyk),
        "v_ldh" => value(v_ldh),
        "v_burn" => value(v_burn),
        "v_lp" => value(v_lp),
    )
    return results_mo
end

end #module

max_stars_repo_path: deprecated/Mos model.jl

max_issues_repo_path: deprecated/Mos model.jl

max_forks_repo_path: deprecated/Mos model.jl


text:
struct Biped2{T}
    n::Int
    m::Int
    d::Int

    l1::T
    l2::T
end

function f(model::Biped2, x, u, w)
    biped_f(x) + biped_g(x) * u
end

n, m, d = 10, 4, 0
model = Biped2(n, m, d, 0.2755, 0.288)

function kinematics(model::Biped2,q)
	θ1 = π - (q[1] + q[5])
	θ2 = π - (q[2] + q[5])

	pe = [model.l1*sin(θ1) + model.l2*sin(θ1 - q[3]),
		  model.l1*cos(θ1) + model.l2*cos(θ1 - q[3])]

	pf = [(model.l1*sin(θ2) + model.l2*sin(θ2 - q[4])),
		  pe[2] - (model.l1*cos(θ2) + model.l2*cos(θ2 - q[4]))]

    return pf
end

function pe(model::Biped2,q)
	θ1 = π - (q[1] + q[5])
	θ2 = π - (q[2] + q[5])

	pe = [-1.0*(model.l1*sin(θ1) + model.l2*sin(θ1 - q[3])),
		  model.l1*cos(θ1) + model.l2*cos(θ1 - q[3])]

    return pe
end

function transformation_to_urdf_left_pinned(model::Biped2, x)
	robot_world_angle = (x[1] + x[3] + x[5]) - pi / 2
	robot_position = [robot_world_angle; x[3]; -x[1] + pi / 2; pi / 2 - x[2]; x[4]]
	return robot_position
end

function transformation_to_urdf_right_pinned(model::Biped2, x)
	robot_world_angle = (x[1] + x[3] + x[5]) - pi / 2
	robot_position = [robot_world_angle; -x[3]; x[1] - pi / 2; x[2] - pi / 2; -x[4]]
	return robot_position
end

function biped_f(in1)
	dq1 = in1[6]
	dq2 = in1[7]
	dq3 = in1[8]
	dq4 = in1[9]
	dq5 = in1[10]
	q1 = in1[1]
	q2 = in1[2]
	q3 = in1[3]
	q4 = in1[4]
	q5 = in1[5]
	t2 = q3.*2.0
	t3 = q4.*2.0
	t4 = q1.*2.0
	t5 = q2.*2.0
	t6 = cos(t2)
	t7 = cos(t3)
	t8 = t2+t4-t5
	t9 = cos(t8)
	t10 = t3-t4+t5
	t11 = cos(t10)
	t12 = t4-t5
	t13 = cos(t12)
	t14 = t2-t3
	t15 = cos(t14)
	t16 = t2+t3
	t17 = cos(t16)
	t18 = dq1.^2
	t19 = q1-q2
	t20 = sin(t19)
	t21 = q1+q3+q5
	t22 = q1+q5
	t23 = sin(t22)
	t24 = dq2.^2
	t25 = q1-q2+q3-q4
	t26 = sin(t25)
	t27 = dq3.^2
	t28 = dq4.^2
	t29 = q1-q2+q3
	t30 = sin(t29)
	t31 = -q1+q2+q4
	t32 = sin(t31)
	t33 = sin(t21)
	t34 = cos(t21)
	t35 = t34.*(5.6e1./3.65e2)
	t36 = cos(t22)
	t37 = t36.*1.372353673723537e-1
	t38 = t35+t37
	t39 = cos(t38)
	t40 = q1+q3
	t41 = sin(t40)
	t42 = sin(q3)
	t43 = sin(q1)
	t44 = sin(q4)
	t45 = q2+q4+q5
	t46 = sin(t45)
	t47 = t18.*t26.*1.35164722e-2
	t48 = t26.*t27.*1.35164722e-2
	t49 = dq5.^2
	t50 = dq1.*dq3.*t26.*2.70329444e-2
	t51 = dq1.*dq5.*t26.*2.70329444e-2
	t52 = dq3.*dq5.*t26.*2.70329444e-2
	t53 = t2-t3+t4-t5
	t54 = cos(t53)
	t55 = cos(t4)
	t56 = cos(t5)
	t57 = t9.*4.646673597987376e90
	t58 = t11.*9.306859101949199e89
	t59 = t2-t3+t4
	t60 = cos(t59)
	t61 = t2+t3+t4
	t62 = cos(t61)
	t63 = -t2+t3+t5
	t64 = cos(t63)
	t65 = t2+t3+t5
	t66 = cos(t65)
	t67 = t2+t4
	t68 = cos(t67)
	t69 = t68.*5.503777366439374e90
	t70 = -t3+t4
	t71 = cos(t70)
	t72 = t71.*5.905084341517107e90
	t73 = t3+t4
	t74 = cos(t73)
	t75 = t74.*5.905084341517107e90
	t76 = -t2+t5
	t77 = cos(t76)
	t78 = t77.*2.948255604070494e91
	t79 = t2+t5
	t80 = cos(t79)
	t81 = t80.*2.948255604070494e91
	t82 = t3+t5
	t83 = cos(t82)
	t84 = t83.*6.340172977765353e89
	t85 = t15.*1.943890255716759e92
	t86 = t17.*1.943890255716759e92
	t89 = t54.*4.996295833675426e88
	t90 = t55.*1.025217124508562e92
	t91 = t56.*5.896511208140987e91
	t92 = t6.*3.89576798740676e93
	t93 = t7.*4.982848297232157e92
	t94 = t60.*3.170086488882676e89
	t95 = t62.*3.170086488882676e89
	t96 = t64.*3.170086488882676e89
	t97 = t66.*3.170086488882676e89
	t98 = t13.*8.655599649991668e91
	t87 = t57+t58+t69+t72+t75+t78+t81+t84+t85+t86-t89-t90-t91-t92-t93-t94-t95-t96-t97-t98+4.920249036048467e93
	t88 = 1.0./t87
	t99 = q1+t2
	t100 = cos(t99)
	t101 = q1+t3
	t102 = cos(t101)
	t103 = -q1+t3+t5
	t104 = cos(t103)
	t105 = t11.*2.113657637960732e66
	t106 = t15.*5.408856160038117e68
	t107 = t17.*5.408856160038117e68
	t108 = q1+t2-t3-t5
	t109 = cos(t108)
	t110 = cos(q1)
	t111 = q1+t2-t5
	t112 = cos(t111)
	t113 = q1+t2-t3
	t114 = cos(t113)
	t115 = q1+t2+t3
	t116 = cos(t115)
	t117 = t33.*2.24837916075e1
	t118 = t23.*1.79439657183e1
	t119 = t23.*t39.*3.062976868808966
	t120 = t26.*t49.*1.35164722e-2
	t121 = t24.*t32.*1.2090221075e-2
	t122 = t28.*t32.*1.2090221075e-2
	t123 = t33.*t39.*3.424308078378082
	t124 = dq2.*dq4.*t32.*2.418044215e-2
	t125 = dq2.*dq5.*t32.*2.418044215e-2
	t126 = dq4.*dq5.*t32.*2.418044215e-2
	t127 = q2+q5
	t128 = sin(t127)
	t129 = t18.*t20.*7.1681372485e-2
	t130 = t20.*t49.*7.1681372485e-2
	t131 = t18.*t30.*8.013743276e-2
	t132 = t27.*t30.*8.013743276e-2
	t133 = t30.*t49.*8.013743276e-2
	t134 = dq1.*dq5.*t20.*1.4336274497e-1
	t135 = dq1.*dq3.*t30.*1.6027486552e-1
	t136 = dq1.*dq5.*t30.*1.6027486552e-1
	t137 = dq3.*dq5.*t30.*1.6027486552e-1
	t138 = cos(t19)
	t139 = q1-t5
	t140 = cos(t139)
	t141 = t100.*1.158914996825368e69
	t142 = q1-t3
	t143 = cos(t142)
	t144 = t102.*7.921782016575439e67
	t145 = q2-t2
	t146 = cos(t145)
	t147 = q2+t2
	t148 = cos(t147)
	t149 = q2+t3
	t150 = cos(t149)
	t151 = -q2+t2-t3+t4
	t152 = cos(t151)
	t153 = -q2+t2+t4
	t154 = cos(t153)
	t155 = t104.*9.035967085704604e65
	t156 = t9.*1.055294491289711e67
	t157 = q1+q2+t2+t3
	t158 = cos(t157)
	t159 = -q2+t4
	t160 = cos(t159)
	t161 = q1+q2+t2
	t162 = cos(t161)
	t163 = q1+q2+t3
	t164 = cos(t163)
	t165 = q1+q2
	t166 = cos(t165)
	t167 = q1-q2+t2-t3
	t168 = cos(t167)
	t169 = cos(q2)
	t170 = q1-q2+t2
	t171 = cos(t170)
	t172 = t112.*8.403663651465465e67
	t173 = -q1+q2+t3
	t174 = cos(t173)
	t175 = q2+t3-t4
	t176 = cos(t175)
	t177 = q2-t2+t3
	t178 = cos(t177)
	t179 = q2+t2+t3
	t180 = cos(t179)
	t181 = t24.*t26.*1.35164722e-2
	t182 = t26.*t28.*1.35164722e-2
	t183 = t24.*t30.*8.013743276e-2
	t184 = t41.*t49.*5.57372816e-2
	t185 = dq2.*dq4.*t26.*2.70329444e-2
	t186 = dq2.*dq5.*t26.*2.70329444e-2
	t187 = dq4.*dq5.*t26.*2.70329444e-2
	t188 = dq2.*dq5.*t30.*1.6027486552e-1
	t189 = cos(q3)
	t190 = cos(q4)
	t191 = t27.*t42.*5.6337833244e-1
	t192 = t43.*t49.*4.98559126e-2
	t193 = t20.*t24.*7.1681372485e-2
	t194 = dq1.*dq3.*t42.*1.12675666488
	t195 = dq3.*dq5.*t42.*1.12675666488
	t196 = dq2.*dq5.*t20.*1.4336274497e-1
	t212 = t32.*t49.*1.2090221075e-2
	t197 = -t117-t118-t119+t120-t121-t122-t123-t124-t125-t126+t130+t133+t181+t182+t183+t184+t185+t186+t187+t188+t191+t192+t193+t194+t195+t196-t212
	t198 = t138.*9.46867650868678e68
	t199 = t143.*7.921782016575439e67
	t200 = t146.*6.531763836093252e68
	t201 = t148.*6.531763836093252e68
	t202 = t150.*1.407002795889341e68
	t203 = t154.*8.243393381131519e66
	t204 = t158.*7.19564337544974e66
	t205 = t166.*6.564490376243981e67
	t206 = t168.*7.110395898963896e67
	t207 = t176.*1.683176403400424e67
	t208 = t189.^2
	t209 = t190.^2
	t210 = t24.*t44.*1.2090221075e-2
	t211 = t44.*t49.*1.2090221075e-2
	t213 = dq2.*dq5.*t44.*2.418044215e-2
	t214 = t110.^2
	t215 = t169.^2
	t250 = t214.*2.050434249017124e92
	t251 = t215.*1.179302241628197e92
	t252 = t208.*7.79153597481352e93
	t253 = t209.*9.965696594464315e92
	t216 = t57+t58+t69+t72+t75+t78+t81+t84+t85+t86-t89-t94-t95-t96-t97-t98-t250-t251-t252-t253+9.475788677710709e93
	t217 = 1.0./t216
	t218 = t146.*1.632940959023313e68
	t219 = t148.*1.632940959023313e68
	t220 = q2-q4
	t221 = cos(t220)
	t222 = t150.*3.517506989723352e67
	t223 = -q2-q4+t2+t4
	t224 = cos(t223)
	t225 = t54.*2.835219413054969e65
	t226 = t154.*2.06084834528288e66
	t227 = q2+q4-t4
	t228 = cos(t227)
	t229 = t228.*2.276934412542559e67
	t230 = q2+q4-t2
	t231 = cos(t230)
	t232 = q2+q4+t2
	t233 = cos(t232)
	t234 = q2+q4
	t235 = cos(t234)
	t236 = t235.*1.876551142359231e68
	t237 = -q2+q4+t2+t4
	t238 = cos(t237)
	t239 = t238.*1.339327292724032e66
	t240 = -q2+q4+t4
	t241 = cos(t240)
	t242 = t176.*4.207941008501061e66
	t243 = q2-q4+t2
	t244 = cos(t243)
	t245 = t244.*9.246973441819591e67
	t246 = -q2+q4+t2
	t247 = cos(t246)
	t248 = t247.*9.246973441819591e67
	t313 = t18.*t42.*5.6337833244e-1
	t314 = t42.*t49.*5.6337833244e-1
	t315 = dq1.*dq5.*t42.*1.12675666488
	t249 = -t117+t120-t123+t133+t181+t182+t183+t184+t185+t186+t187+t188-t313-t314-t315
	t254 = q1-q3
	t255 = cos(t254)
	t256 = q1+q3-t3-t5
	t257 = cos(t256)
	t258 = q1+q3-t5
	t259 = cos(t258)
	t260 = q1+q3-t3
	t261 = cos(t260)
	t262 = q1+q3+t3
	t263 = cos(t262)
	t264 = cos(t40)
	t265 = -q1+q3+t3+t5
	t266 = cos(t265)
	t267 = -q1+q3+t5
	t268 = cos(t267)
	t269 = q1-q3+t3
	t270 = cos(t269)
	t271 = -q1+q3+t3
	t272 = cos(t271)
	t273 = t54.*1.020618132398732e66
	t274 = sin(q5)
	t275 = t18.*t43.*4.98559126e-2
	t276 = t18.*t41.*5.57372816e-2
	t277 = t27.*t41.*5.57372816e-2
	t278 = dq1.*dq5.*t43.*9.97118252e-2
	t279 = dq1.*dq3.*t41.*1.114745632e-1
	t280 = dq1.*dq5.*t41.*1.114745632e-1
	t281 = dq3.*dq5.*t41.*1.114745632e-1
	t282 = t100.*2.89728749206342e68
	t283 = q3+t4
	t284 = cos(t283)
	t285 = t143.*1.98044550414386e67
	t286 = t102.*1.98044550414386e67
	t287 = q3-t3
	t288 = cos(t287)
	t289 = q3+t3
	t290 = cos(t289)
	t291 = t104.*2.258991771426151e65
	t292 = t60.*1.798910843862435e66
	t293 = t62.*1.798910843862435e66
	t294 = t259.*1.879228357125597e67
	t295 = t264.*2.556485952983844e68
	t296 = t266.*2.020624133207483e65
	t297 = q3-t3+t4-t5
	t298 = cos(t297)
	t299 = t298.*2.536048533430662e65
	t300 = t189.*2.973264907774096e68
	t301 = t112.*2.100915912866366e67
	t302 = q3+t4-t5
	t303 = cos(t302)
	t304 = t270.*1.395077797426718e67
	t305 = t272.*1.395077797426718e67
	t306 = q3-t3+t4
	t307 = cos(t306)
	t308 = t307.*1.609090706117211e66
	t309 = q3+t3+t4
	t310 = cos(t309)
	t311 = t310.*1.609090706117211e66
	t312 = t7.*-4.083328958509943e67-t9.*2.63682128314895e67-t68.*3.123197055163642e67-t109.*2.258991771426151e65-t110.*3.627853379244088e68-t114.*1.559651403293754e67-t116.*1.559651403293754e67-t140.*2.100915912866366e67+t225-t255.*2.591567220985299e68-t257.*2.020624133207483e65-t261.*1.374871556094643e67-t263.*1.374871556094643e67-t268.*1.879228357125597e67+t282-t284.*2.793638924342803e67+t285+t286-t288.*1.590027654354345e67-t290.*1.590027654354345e67+t291+t292+t293+t294+t295+t296+t299+t300+t301-t303.*2.358585271128362e67+t304+t305+t308+t311+3.782473003615523e68
	t368 = t46.*4.305084165e-1
	t369 = t128.*2.5524292707
	t370 = t28.*t44.*1.2090221075e-2
	t371 = t18.*t32.*1.2090221075e-2
	t372 = dq2.*dq4.*t44.*2.418044215e-2
	t373 = dq4.*dq5.*t44.*2.418044215e-2
	t374 = dq1.*dq5.*t32.*2.418044215e-2
	t316 = t47+t48+t50+t51+t52+t120+t129+t130+t131+t132+t133+t134+t135+t136+t137-t212-t368-t369-t370-t371-t372-t373-t374
	t317 = t138.*2.958961408964619e67
	t318 = t100.*3.621609365079275e67
	t319 = t143.*2.475556880179825e66
	t320 = t102.*2.475556880179825e66
	t321 = q1+q2+q3
	t322 = cos(t321)
	t323 = t104.*2.823739714282689e64
	t324 = q1+q2-q3+t3
	t325 = cos(t324)
	t326 = t325.*2.011363382646513e65
	t327 = q1-q2+q3-t3
	t328 = cos(t327)
	t329 = -q1+q2+q3+t3
	t330 = cos(t329)
	t331 = t330.*1.987534567942932e66
	t332 = t257.*2.525780166509353e64
	t333 = q1+q2-q3
	t334 = cos(t333)
	t335 = cos(t29)
	t336 = t335.*2.079104639561109e67
	t337 = -q1+q2+q3
	t338 = cos(t337)
	t339 = t162.*2.051403242576244e66
	t340 = t164.*2.248638554828044e65
	t341 = t261.*1.718589445118304e66
	t342 = t263.*1.718589445118304e66
	t343 = t109.*2.823739714282689e64
	t344 = t110.*4.53481672405511e67
	t345 = t112.*2.626144891082958e66
	t346 = t174.*2.882162479711211e66
	t347 = t114.*1.949564254117192e66
	t348 = t116.*1.949564254117192e66
	t349 = q1+q2+q3+t3
	t350 = cos(t349)
	t351 = t140.*2.626144891082958e66
	t352 = t158.*2.248638554828044e65
	t353 = t166.*2.051403242576244e66
	t354 = t168.*2.221998718426218e66
	t355 = t171.*2.356702732209462e67
	t356 = t255.*3.239459026231624e67
	t357 = t259.*2.349035446406996e66
	t358 = t264.*3.195607441229805e67
	t359 = t266.*2.525780166509353e64
	t360 = t268.*2.349035446406996e66
	t361 = t270.*1.743847246783397e66
	t362 = t272.*1.743847246783397e66
	t363 = t328.*1.955833961275048e66
	t364 = t334.*1.834940238083621e66
	t365 = t338.*2.108024684167879e67
	t366 = t350.*2.011363382646513e65
	t367 = t317+t318+t319+t320-t322.*1.834940238083621e66+t323-t326-t331-t332-t336-t339-t340-t341-t342-t343-t344+t345-t346-t347-t348-t351+t352+t353+t354-t355-t356+t357+t358+t359-t360+t361+t362+t363+t364+t365+t366
	t375 = q1+q2+q4
	t376 = cos(t375)
	t377 = q1+q2-q4+t2
	t378 = cos(t377)
	t379 = q1-q2+q4+t2
	t380 = cos(t379)
	t381 = q1+q2-q4
	t382 = cos(t381)
	t383 = q1-q2+q4
	t384 = cos(t383)
	t385 = cos(t31)
	t386 = q1-q2-q4+t2
	t387 = cos(t386)
	t388 = q1+q2+q4+t2
	t389 = cos(t388)
	t390 = q4+t4
	t391 = cos(t390)
	t392 = q4+t5
	t393 = cos(t392)
	t394 = q4+t2
	t395 = cos(t394)
	t396 = -q4+t2+t4-t5
	t397 = cos(t396)
	t398 = t397.*1.680965282298581e66
	t399 = -q4+t2+t4
	t400 = cos(t399)
	t401 = t400.*9.315304957758298e65
	t402 = t64.*1.798910843862435e66
	t403 = t66.*1.798910843862435e66
	t404 = t68.*1.590481329583487e66
	t405 = -q4+t4
	t406 = cos(t405)
	t407 = -q4+t2
	t408 = cos(t407)
	t409 = t190.*1.704625264978413e69
	t410 = q4-t4+t5
	t411 = cos(t410)
	t412 = q4+t2+t4
	t413 = cos(t412)
	t414 = t413.*9.315304957758298e65
	t415 = q4-t2+t5
	t416 = cos(t415)
	t417 = t416.*1.0665512025487e67
	t418 = q4+t2+t5
	t419 = cos(t418)
	t420 = t419.*1.0665512025487e67
	t421 = t47+t48+t50+t51+t52+t120+t210+t211-t212+t213-t368-t371-t374
	t422 = t54.*3.990216591296235e80
	t423 = t138.*4.16436796838023e82
	t424 = q1+q2-q3-q4
	t425 = cos(t424)
	t426 = t425.*1.67831086311752e81
	t427 = cos(t25)
	t428 = t427.*1.462987332903369e82
	t429 = q1-q2-q3+q4
	t430 = cos(t429)
	t431 = t430.*1.658427753522667e82
	t432 = t328.*2.752591593499022e81
	t433 = t158.*3.164677424254534e80
	t434 = t335.*2.9260796499974e82
	t435 = t162.*2.887093399641563e81
	t436 = t385.*2.165299823858942e82
	t437 = t164.*3.164677424254534e80
	t438 = q1+q2+q3+q4
	t439 = cos(t438)
	t440 = t439.*1.531726456607094e81
	t441 = t166.*2.887093399641563e81
	t442 = t387.*1.662561080786584e82
	t443 = t168.*3.127185187600635e81
	t444 = t171.*3.316764233312988e82
	t445 = t174.*4.056283084265199e81
	t446 = q1+q2+q3-q4
	t447 = cos(t446)
	t448 = t447.*1.67831086311752e81
	t449 = q1+q2-q3+q4
	t450 = cos(t449)
	t451 = t450.*1.531726456607094e81
	t452 = q1-q2+q3+q4
	t453 = cos(t452)
	t454 = t453.*1.631976306212276e82
	t455 = -q1+q2+q3+q4
	t456 = cos(t455)
	t457 = t456.*1.487128499210077e82
	t458 = t322.*2.582448803900164e81
	t459 = t325.*2.830742306435468e80
	t460 = t330.*2.797206231116802e81
	t461 = t334.*2.582448803900164e81
	t462 = t338.*2.966780994408209e82
	t463 = t350.*2.830742306435468e80
	t464 = t376.*1.71242014023588e81
	t465 = t378.*1.876296718113235e81
	t466 = t380.*1.854068051125159e82
	t467 = t382.*1.876296718113235e81
	t468 = t384.*2.404918296068613e82
	t469 = t389.*1.71242014023588e81
	t470 = t423+t426+t428+t431+t432+t433-t434-t435-t436-t437+t440+t441+t442+t443-t444-t445-t448-t451-t454-t457-t458-t459-t460+t461+t462+t463-t464-t465-t466+t467+t468+t469
	t492 = t274.*1.775268612
	t471 = t47+t48+t50+t51+t52+t117+t118+t119+t121+t122+t123+t124+t125+t126+t129+t131+t132+t134+t135+t136+t137-t181-t182-t183-t185-t186-t187-t188-t191-t193-t194-t195-t196+t275+t276+t277+t278+t279+t280+t281-t368-t369-t370-t371-t372-t373-t374-t492
	t472 = t138.*2.367169127171695e68
	t473 = t221.*2.085484849158977e68
	t474 = t224.*1.222349859852932e66
	t475 = t152.*2.258991771426151e65
	t476 = t158.*1.798910843862435e66
	t477 = t160.*3.838848983032198e67
	t478 = t162.*1.641122594060995e67
	t479 = t385.*1.230830448468477e68
	t480 = t164.*1.798910843862435e66
	t481 = t166.*1.641122594060995e67
	t482 = t387.*9.450565589682945e67
	t483 = t168.*1.777598974740974e67
	t484 = t169.*3.629158332897017e68
	t485 = t171.*1.88536218576757e68
	t486 = t174.*2.305729983768969e67
	t487 = t241.*2.494834337222995e67
	t488 = t178.*1.559651403293754e67
	t489 = t180.*1.559651403293754e67
	t490 = t389.*9.733981529711172e66
	t491 = t218+t219+t222+t226+t229-t231.*8.305405998517524e67-t233.*8.305405998517524e67+t236+t239+t242+t245+t248-t376.*9.733981529711172e66-t378.*1.0665512025487e67-t380.*1.053915668265494e68+t382.*1.0665512025487e67+t384.*1.367037780294489e68+t472-t473-t474-t475+t476-t477-t478-t479-t480+t481+t482+t483-t484-t485-t486-t487-t488-t489+t490
	t493 = t100.*4.002547854592237e61
	t494 = t143.*2.735947994619766e60
	t495 = t102.*2.735947994619766e60
	t496 = t104.*3.120754392869666e58
	t497 = t257.*2.791454010505172e58
	t498 = t261.*1.899359042642789e60
	t499 = t263.*1.899359042642789e60
	t500 = t109.*3.120754392869666e58
	t501 = t110.*5.011810805674362e61
	t502 = t112.*2.902375585010773e60
	t503 = t114.*2.154628905576533e60
	t504 = t116.*2.154628905576533e60
	t505 = t140.*2.902375585010773e60+t255.*3.580201084221302e61-t259.*2.596118421008013e60-t264.*3.531736976202944e61-t266.*2.791454010505172e58+t268.*2.596118421008013e60-t270.*1.92727358274784e60-t272.*1.92727358274784e60-t493-t494-t495-t496+t497+t498+t499+t500+t501-t502+t503+t504
	t506 = t146.*2.819846300254332e60
	t507 = t148.*2.819846300254332e60
	t508 = t150.*6.074211695334503e59
	t509 = t224.*2.110816506179317e58
	t510 = t152.*3.900942991087083e57
	t511 = t154.*3.558778748073662e58
	t512 = t160.*6.62912155042827e59
	t513 = t231.*1.434220156437109e60
	t514 = t233.*1.434220156437109e60
	t515 = t169.*6.267016968070755e60
	t516 = t176.*7.266488613040416e58
	t517 = t178.*2.693286131970666e59
	t518 = t180.*2.693286131970666e59
	vec([dq1 dq2 dq3 dq4 dq5 t88.*t491.*(t46.*(-4.305084165e-1)+t47+t48+t50+t51+t52+t120+t210+t211+t213-t18.*t32.*1.2090221075e-2-t32.*t49.*1.2090221075e-2-dq1.*dq5.*t32.*2.418044215e-2).*5.62949953421312e25+t88.*t197.*(t6.*-1.069433799329253e70-t7.*1.519184392954739e69-t9.*9.491990641306088e67-t13.*1.96575172256594e68+t60.*7.19564337544974e66+t62.*7.19564337544974e66-t68.*1.249278822065457e68+t100.*2.317829993650736e69+t102.*1.584356403315088e68+t104.*1.807193417140921e66+t105+t106+t107-t109.*1.807193417140921e66-t110.*2.90228270339527e69+t112.*1.680732730293093e68-t114.*1.247721122635003e68-t116.*1.247721122635003e68-t140.*1.680732730293093e68+t143.*1.584356403315088e68+t273+1.475433566123788e70).*1.40737488355328e25-t88.*(t46.*(-4.305084165e-1)+t47+t48+t50+t51+t52+t120-t128.*2.5524292707+t129+t130+t131+t132+t133+t134+t135+t136+t137-t18.*t32.*1.2090221075e-2-t28.*t44.*1.2090221075e-2-t32.*t49.*1.2090221075e-2-dq1.*dq5.*t32.*2.418044215e-2-dq2.*dq4.*t44.*2.418044215e-2-dq4.*dq5.*t44.*2.418044215e-2).*(t6.*-1.069433799329253e70-t7.*1.355851234614342e69-t13.*1.96575172256594e68-t54.*1.13469632823255e65+t105+t106+t107-t109.*9.035967085704604e65-t110.*1.451141351697635e69-t114.*6.238605613175014e67-t116.*6.238605613175014e67-t140.*8.403663651465465e67+t141+t144-t152.*9.035967085704604e65+t155+t156-t160.*1.535539593212879e68-t162.*6.564490376243981e67-t164.*7.19564337544974e66-t169.*1.451663333158807e69-t171.*7.541448743070279e68+t172-t174.*9.222919935075876e67-t178.*6.238605613175014e67-t180.*6.238605613175014e67+t198+t199+t200+t201+t202+t203+t204+t205+t206+t207+1.324134645979167e70).*1.40737488355328e25-t88.*t249.*t312.*5.62949953421312e25+t88.*(t6.*-1.069433799329253e70-t7.*1.355851234614342e69-t13.*1.96575172256594e68+t105+t106+t107-t109.*9.035967085704604e65-t110.*1.451141351697635e69-t114.*6.238605613175014e67-t116.*6.238605613175014e67+t141+t144+t155+t156+t172-cos(q1-q2.*2.0).*8.403663651465465e67+cos(q1-q4.*2.0).*7.921782016575439e67-cos(q2.*-2.0+t2-t3+t4).*1.13469632823255e65+1.324134645979167e70).*(t46.*(-4.305084165e-1)+t47+t48+t50+t51+t52+t117+t118+t119+t121+t122+t123+t124+t125+t126-t128.*2.5524292707+t129+t131+t132+t134+t135+t136+t137-t274.*1.775268612+t275+t276+t277+t278+t279+t280+t281-t20.*t24.*7.1681372485e-2-t18.*t32.*1.2090221075e-2-t24.*t26.*1.35164722e-2-t24.*t30.*8.013743276e-2-t26.*t28.*1.35164722e-2-t27.*t42.*5.6337833244e-1-t28.*t44.*1.2090221075e-2-dq2.*dq5.*t20.*1.4336274497e-1-dq2.*dq4.*t26.*2.70329444e-2-dq2.*dq5.*t26.*2.70329444e-2-dq4.*dq5.*t26.*2.70329444e-2-dq2.*dq5.*t30.*1.6027486552e-1-dq1.*dq5.*t32.*2.418044215e-2-dq1.*dq3.*t42.*1.12675666488-dq2.*dq4.*t44.*2.418044215e-2-dq3.*dq5.*t42.*1.12675666488-dq4.*dq5.*t44.*2.418044215e-2).*1.40737488355328e25 t217.*(t46.*(-4.305084165e-1)+t47+t48+t50+t51+t52+t120+t210+t211-t212+t213-t18.*t32.*1.2090221075e-2-dq1.*dq5.*t32.*2.418044215e-2).*(t11.*-5.28131009027995e66-t83.*3.59782168772487e66-t152.*2.258991771426151e65-t160.*3.838848983032198e67-t169.*3.629158332897017e68-t178.*1.559651403293754e67-t180.*1.559651403293754e67-t208.*2.651967948720325e69-t214.*5.925343947388025e67+t218+t219-t221.*2.085484849158977e68+t222-t224.*1.222349859852932e66+t225+t226+t229-t231.*8.305405998517524e67-t233.*8.305405998517524e67+t236+t239-t241.*2.494834337222995e67+t242+t245+t248-t391.*1.73521011604649e67-t393.*2.1331024050974e67-t395.*6.766040747477027e68+t398+t401+t402+t403+t404-t406.*1.73521011604649e67-t408.*6.766040747477027e68+t409-t411.*3.131221120289946e67+t414+t417+t420+3.001895035993504e69).*5.62949953421312e25+t217.*(t13.*-1.96575172256594e68-t54.*1.13469632823255e65+t105+t106+t107-t152.*9.035967085704604e65+t156-t160.*1.535539593212879e68-t169.*1.451663333158807e69-t178.*6.238605613175014e67-t180.*6.238605613175014e67+t200+t201+t202+t203+t207-t208.*2.138867598658505e70-t209.*2.711702469228684e69+2.529153568769854e70).*(t46.*(-4.305084165e-1)+t47+t48+t50+t51+t52+t117+t118+t119+t121+t122+t123+t124+t125+t126-t128.*2.5524292707+t129+t131+t132+t134+t135+t136+t137-t181-t182-t183-t185-t186-t187-t188-t191-t193-t194-t195-t196-t274.*1.775268612+t275+t276+t277+t278+t279+t280+t281-t18.*t32.*1.2090221075e-2-t28.*t44.*1.2090221075e-2-dq1.*dq5.*t32.*2.418044215e-2-dq2.*dq4.*t44.*2.418044215e-2-dq4.*dq5.*t44.*2.418044215e-2).*1.40737488355328e25-t217.*t249.*t367.*4.503599627370496e26-t217.*t316.*(t11.*-1.901158272315907e67-t13.*1.96575172256594e68+t64.*7.19564337544974e66+t66.*7.19564337544974e66+t68.*6.36192531833395e66-t83.*1.439128675089948e67+t106+t107+t146.*1.30635276721865e69+t148.*1.30635276721865e69+t150.*2.814005591778681e68-t152.*1.807193417140921e66+t154.*1.648678676226304e67+t156-t160.*3.071079186425758e68-t169.*2.903326666317613e69+t176.*3.366352806800849e67-t178.*1.247721122635003e68-t180.*1.247721122635003e68-t208.*3.199654778146635e70-t209.*2.711702469228684e69-t214.*2.37013757895521e68+t273+3.729911583167255e70).*1.40737488355328e25+t197.*t217.*(t13.*-1.96575172256594e68-t54.*1.13469632823255e65+t105+t106+t107-t109.*9.035967085704604e65-t110.*1.451141351697635e69-t114.*6.238605613175014e67-t116.*6.238605613175014e67-t140.*8.403663651465465e67+t141+t144-t152.*9.035967085704604e65+t155+t156-t160.*1.535539593212879e68-t162.*6.564490376243981e67-t164.*7.19564337544974e66-t169.*1.451663333158807e69-t171.*7.541448743070279e68+t172-t174.*9.222919935075876e67-t178.*6.238605613175014e67-t180.*6.238605613175014e67+t198+t199+t200+t201+t202+t203+t204+t205+t206+t207-t208.*2.138867598658505e70-t209.*2.711702469228684e69+2.529153568769854e70).*1.40737488355328e25 t88.*t197.*t312.*-5.62949953421312e25+t88.*t316.*t367.*4.503599627370496e26-t88.*t421.*t470.*3.2e11+t88.*t471.*t505.*4.074974036435136e32+t88.*t249.*(t7.*-9.718123450432042e82-t9.*3.710996046322563e82+t11.*3.192553938622998e80-t13.*2.969150865079421e82-t55.*3.51682660384109e82+t60.*2.531741939403627e81+t62.*2.531741939403627e81-t68.*4.395509091824878e82+t71.*2.025635079020705e81+t74.*2.025635079020705e81+t189.*8.368996706703245e83-t284.*7.863394511673726e82-t288.*4.475529969786883e82-t290.*4.475529969786883e82+t298.*7.138342018844889e80-t303.*6.638827342609521e82+t307.*4.529187690296749e81+t310.*4.529187690296749e81+t422+9.038595841607977e83).*4.0e10 t88.*t421.*(t6.*-1.257359213783462e85+t9.*1.402621882323949e82-t11.*7.432783173316508e81-t13.*2.612736448493452e83-t55.*3.186269261421593e83-t56.*1.779891674226747e83+t64.*2.531741939403627e81+t66.*2.531741939403627e81+t68.*1.710517335808237e82+t77.*8.899458371133735e82+t80.*8.899458371133735e82-t83.*5.063483878807254e81+t190.*4.798093567601945e84-t391.*4.884182270022806e82-t393.*6.004149497962351e82-t395.*1.904471161819446e84+t397.*4.731496636864144e81+t400.*2.622025246037675e81-t406.*4.884182270022806e82-t408.*1.904471161819446e84-t411.*8.813603919095266e82+t413.*2.622025246037675e81+t416.*3.002074748981175e82+t419.*3.002074748981175e82+t422+1.581779989358623e85).*-4.0e10+t88.*t316.*(t6.*-1.325983974360163e69-t11.*5.28131009027995e66-t55.*2.962671973694012e67-t83.*3.59782168772487e66-t152.*2.258991771426151e65-t160.*3.838848983032198e67-t169.*3.629158332897017e68-t178.*1.559651403293754e67-t180.*1.559651403293754e67+t218+t219-t221.*2.085484849158977e68+t222-t224.*1.222349859852932e66+t225+t226+t229-t231.*8.305405998517524e67-t233.*8.305405998517524e67+t236+t239-t241.*2.494834337222995e67+t242+t245+t248-t391.*1.73521011604649e67-t393.*2.1331024050974e67-t395.*6.766040747477027e68+t398+t401+t402+t403+t404-t406.*1.73521011604649e67-t408.*6.766040747477027e68+t409-t411.*3.131221120289946e67+t414+t417+t420+1.646284341896401e69).*5.62949953421312e25+t88.*t471.*(t221.*3.601322328061862e60-t228.*3.93192726512915e59-t235.*3.240524874325339e60-t238.*2.312819144102071e58+t241.*4.308208044321151e59-t244.*1.596814857535384e60-t247.*1.596814857535384e60-t506-t507-t508+t509+t510-t511+t512+t513+t514+t515-t516+t517+t518).*3.259979229148109e33-t88.*t197.*t491.*5.62949953421312e25+t88.*t249.*t470.*3.2e11 t88.*t197.*(t6.*-1.069433799329253e70-t7.*1.355851234614342e69-t13.*1.96575172256594e68-t54.*1.13469632823255e65+t105+t106+t107-t109.*9.035967085704604e65-t110.*1.451141351697635e69-t114.*6.238605613175014e67-t116.*6.238605613175014e67-t140.*8.403663651465465e67+t141+t144+t155+t156+t172+t199+1.324134645979167e70).*-1.40737488355328e25+t88.*t316.*(t6.*-1.069433799329253e70-t7.*1.355851234614342e69-t13.*1.96575172256594e68-t54.*1.13469632823255e65+t105+t106+t107-t152.*9.035967085704604e65+t156-t160.*1.535539593212879e68-t169.*1.451663333158807e69-t178.*6.238605613175014e67-t180.*6.238605613175014e67+t200+t201+t202+t203+t207+1.324134645979167e70).*1.40737488355328e25+t88.*t421.*(t221.*3.601322328061862e60-t228.*3.93192726512915e59-t235.*3.240524874325339e60-t238.*2.312819144102071e58+t241.*4.308208044321151e59-t244.*1.596814857535384e60-t247.*1.596814857535384e60-t506-t507-t508+t509+t510-t511+t512+t513+t514+t515-t516+t517+t518).*3.259979229148109e33-t88.*t249.*t505.*4.074974036435136e32-t88.*t471.*(t6.*-1.069433799329253e70-t7.*1.355851234614342e69-t13.*1.96575172256594e68-t54.*1.13469632823255e65+t105+t106+t107+t156+1.324134645979167e70).*1.40737488355328e25])
end

function biped_g(in1)
	q1 = in1[1]
	q2 = in1[2]
	q3 = in1[3]
	q4 = in1[4]
	t2 = q1.*2.0
	t3 = q3.*2.0
	t4 = q2.*2.0
	t5 = q4.*2.0
	t6 = cos(t3)
	t7 = cos(t5)
	t8 = t2+t3-t4
	t9 = cos(t8)
	t10 = -t2+t4+t5
	t11 = cos(t10)
	t12 = t2+t3-t5
	t13 = cos(t12)
	t14 = t2+t3+t5
	t15 = cos(t14)
	t16 = t2-t4
	t17 = cos(t16)
	t18 = t2+t3
	t19 = cos(t18)
	t20 = t3-t5
	t21 = cos(t20)
	t22 = t3+t5
	t23 = cos(t22)
	t24 = t2+t3-t4-t5
	t25 = cos(t24)
	t26 = cos(t2)
	t27 = cos(t4)
	t28 = t9.*4.646673597987376e90
	t29 = t11.*9.306859101949199e89
	t30 = -t3+t4+t5
	t31 = cos(t30)
	t32 = t3+t4+t5
	t33 = cos(t32)
	t34 = t19.*5.503777366439374e90
	t35 = t2-t5
	t36 = cos(t35)
	t37 = t36.*5.905084341517107e90
	t38 = t2+t5
	t39 = cos(t38)
	t40 = t39.*5.905084341517107e90
	t41 = -t3+t4
	t42 = cos(t41)
	t43 = t42.*2.948255604070494e91
	t44 = t3+t4
	t45 = cos(t44)
	t46 = t45.*2.948255604070494e91
	t47 = t4+t5
	t48 = cos(t47)
	t49 = t48.*6.340172977765353e89
	t50 = t21.*1.943890255716759e92
	t51 = t23.*1.943890255716759e92
	t52 = q1-t4
	t53 = cos(t52)
	t54 = q1+t3
	t55 = cos(t54)
	t56 = q1-t5
	t57 = cos(t56)
	t58 = q1+t5
	t59 = cos(t58)
	t60 = -q1+t4+t5
	t61 = cos(t60)
	t62 = t11.*2.113657637960732e66
	t63 = t21.*5.408856160038117e68
	t64 = t23.*5.408856160038117e68
	t65 = q1+t3-t4-t5
	t66 = cos(t65)
	t67 = cos(q1)
	t68 = q1+t3-t4
	t69 = cos(t68)
	t70 = q1+t3-t5
	t71 = cos(t70)
	t72 = q1+t3+t5
	t73 = cos(t72)
	t76 = t25.*4.996295833675426e88
	t77 = t26.*1.025217124508562e92
	t78 = t27.*5.896511208140987e91
	t79 = t6.*3.89576798740676e93
	t80 = t7.*4.982848297232157e92
	t81 = t13.*3.170086488882676e89
	t82 = t15.*3.170086488882676e89
	t83 = t31.*3.170086488882676e89
	t84 = t33.*3.170086488882676e89
	t85 = t17.*8.655599649991668e91
	t74 = t28+t29+t34+t37+t40+t43+t46+t49+t50+t51-t76-t77-t78-t79-t80-t81-t82-t83-t84-t85+4.920249036048467e93
	t75 = 1.0./t74
	t86 = q1-q2
	t87 = cos(t86)
	t88 = q2-t3
	t89 = cos(t88)
	t90 = q2+t3
	t91 = cos(t90)
	t92 = q2+t5
	t93 = cos(t92)
	t94 = -q2+t2+t3-t5
	t95 = cos(t94)
	t96 = -q2+t2+t3
	t97 = cos(t96)
	t98 = q1+q2+t3+t5
	t99 = cos(t98)
	t100 = -q2+t2
	t101 = cos(t100)
	t102 = q1+q2+t3
	t103 = cos(t102)
	t104 = q1+q2+t5
	t105 = cos(t104)
	t106 = q1+q2
	t107 = cos(t106)
	t108 = q1-q2+t3-t5
	t109 = cos(t108)
	t110 = cos(q2)
	t111 = q1-q2+t3
	t112 = cos(t111)
	t113 = -q1+q2+t5
	t114 = cos(t113)
	t115 = q2-t2+t5
	t116 = cos(t115)
	t117 = q2-t3+t5
	t118 = cos(t117)
	t119 = q2+t3+t5
	t120 = cos(t119)
	t121 = t87.*9.46867650868678e68
	t122 = t55.*1.158914996825368e69
	t123 = t57.*7.921782016575439e67
	t124 = t59.*7.921782016575439e67
	t125 = t89.*6.531763836093252e68
	t126 = t91.*6.531763836093252e68
	t127 = t93.*1.407002795889341e68
	t128 = t97.*8.243393381131519e66
	t129 = t61.*9.035967085704604e65
	t130 = t9.*1.055294491289711e67
	t131 = t99.*7.19564337544974e66
	t132 = t107.*6.564490376243981e67
	t133 = t109.*7.110395898963896e67
	t134 = t69.*8.403663651465465e67
	t135 = t116.*1.683176403400424e67
	t139 = t7.*1.355851234614342e69
	t140 = t17.*1.96575172256594e68
	t382 = t53.*8.403663651465465e67
	t383 = t25.*1.13469632823255e65
	t384 = t6.*1.069433799329253e70
	t385 = t66.*9.035967085704604e65
	t386 = t67.*1.451141351697635e69
	t387 = t71.*6.238605613175014e67
	t388 = t73.*6.238605613175014e67
	t389 = t95.*9.035967085704604e65
	t390 = t101.*1.535539593212879e68
	t391 = t110.*1.451663333158807e69
	t392 = t118.*6.238605613175014e67
	t393 = t120.*6.238605613175014e67
	t136 = t62+t63+t64-t103.*6.564490376243981e67-t105.*7.19564337544974e66-t112.*7.541448743070279e68-t114.*9.222919935075876e67+t121+t122+t123+t124+t125+t126+t127+t128+t129+t130+t131+t132+t133+t134+t135-t139-t140-t382-t383-t384-t385-t386-t387-t388-t389-t390-t391-t392-t393+1.324134645979167e70
	t137 = t75.*t136.*1.40737488355328e25
	t138 = t25.*1.020618132398732e66
	t141 = q1-q3
	t142 = cos(t141)
	t143 = q1+q3-t4-t5
	t144 = cos(t143)
	t145 = q1+q3-t4
	t146 = cos(t145)
	t147 = q1+q3-t5
	t148 = cos(t147)
	t149 = q1+q3+t5
	t150 = cos(t149)
	t151 = q1+q3
	t152 = cos(t151)
	t153 = -q1+q3+t4+t5
	t154 = cos(t153)
	t155 = -q1+q3+t4
	t156 = cos(t155)
	t157 = q1-q3+t5
	t158 = cos(t157)
	t159 = -q1+q3+t5
	t160 = cos(t159)
	t161 = t89.*1.632940959023313e68
	t162 = t91.*1.632940959023313e68
	t163 = q2-q4
	t164 = cos(t163)
	t165 = t93.*3.517506989723352e67
	t166 = -q2-q4+t2+t3
	t167 = cos(t166)
	t168 = t25.*2.835219413054969e65
	t169 = t97.*2.06084834528288e66
	t170 = q2+q4-t2
	t171 = cos(t170)
	t172 = t171.*2.276934412542559e67
	t173 = q2+q4-t3
	t174 = cos(t173)
	t175 = q2+q4+t3
	t176 = cos(t175)
	t177 = q2+q4
	t178 = cos(t177)
	t179 = t178.*1.876551142359231e68
	t180 = -q2+q4+t2+t3
	t181 = cos(t180)
	t182 = t181.*1.339327292724032e66
	t183 = -q2+q4+t2
	t184 = cos(t183)
	t185 = t116.*4.207941008501061e66
	t186 = q2-q4+t3
	t187 = cos(t186)
	t188 = t187.*9.246973441819591e67
	t189 = -q2+q4+t3
	t190 = cos(t189)
	t191 = t190.*9.246973441819591e67
	t192 = t55.*2.89728749206342e68
	t193 = q3+t2
	t194 = cos(t193)
	t195 = t57.*1.98044550414386e67
	t196 = t59.*1.98044550414386e67
	t197 = q3-t5
	t198 = cos(t197)
	t199 = q3+t5
	t200 = cos(t199)
	t201 = t61.*2.258991771426151e65
	t202 = t13.*1.798910843862435e66
	t203 = t15.*1.798910843862435e66
	t204 = t146.*1.879228357125597e67
	t205 = t152.*2.556485952983844e68
	t206 = t154.*2.020624133207483e65
	t207 = q3+t2-t4-t5
	t208 = cos(t207)
	t209 = t208.*2.536048533430662e65
	t210 = cos(q3)
	t211 = t210.*2.973264907774096e68
	t212 = t69.*2.100915912866366e67
	t213 = q3+t2-t4
	t214 = cos(t213)
	t215 = t158.*1.395077797426718e67
	t216 = t160.*1.395077797426718e67
	t217 = q3+t2-t5
	t218 = cos(t217)
	t219 = t218.*1.609090706117211e66
	t220 = q3+t2+t5
	t221 = cos(t220)
	t222 = t221.*1.609090706117211e66
	t223 = t7.*-4.083328958509943e67-t9.*2.63682128314895e67-t19.*3.123197055163642e67-t53.*2.100915912866366e67-t66.*2.258991771426151e65-t67.*3.627853379244088e68-t71.*1.559651403293754e67-t73.*1.559651403293754e67-t142.*2.591567220985299e68-t144.*2.020624133207483e65-t148.*1.374871556094643e67-t150.*1.374871556094643e67-t156.*1.879228357125597e67+t168+t192-t194.*2.793638924342803e67+t195+t196-t198.*1.590027654354345e67-t200.*1.590027654354345e67+t201+t202+t203+t204+t205+t206+t209+t211+t212-t214.*2.358585271128362e67+t215+t216+t219+t222+3.782473003615523e68
	t224 = t53.*2.626144891082958e66
	t225 = t55.*3.621609365079275e67
	t226 = t57.*2.475556880179825e66
	t227 = t59.*2.475556880179825e66
	t228 = q1+q2+q3
	t229 = cos(t228)
	t230 = t61.*2.823739714282689e64
	t231 = q1+q2-q3+t5
	t232 = cos(t231)
	t233 = t232.*2.011363382646513e65
	t234 = q1-q2+q3-t5
	t235 = cos(t234)
	t236 = -q1+q2+q3+t5
	t237 = cos(t236)
	t238 = t237.*1.987534567942932e66
	t239 = t144.*2.525780166509353e64
	t240 = q1+q2-q3
	t241 = cos(t240)
	t242 = q1-q2+q3
	t243 = cos(t242)
	t244 = t243.*2.079104639561109e67
	t245 = -q1+q2+q3
	t246 = cos(t245)
	t247 = t103.*2.051403242576244e66
	t248 = t105.*2.248638554828044e65
	t249 = t148.*1.718589445118304e66
	t250 = t150.*1.718589445118304e66
	t251 = t66.*2.823739714282689e64
	t252 = t67.*4.53481672405511e67
	t253 = t69.*2.626144891082958e66
	t254 = t114.*2.882162479711211e66
	t255 = t71.*1.949564254117192e66
	t256 = t73.*1.949564254117192e66
	t257 = q1+q2+q3+t5
	t258 = cos(t257)
	t259 = t258.*2.011363382646513e65
	t260 = t87.*2.958961408964619e67+t99.*2.248638554828044e65+t107.*2.051403242576244e66+t109.*2.221998718426218e66-t112.*2.356702732209462e67-t142.*3.239459026231624e67+t146.*2.349035446406996e66+t152.*3.195607441229805e67+t154.*2.525780166509353e64-t156.*2.349035446406996e66+t158.*1.743847246783397e66+t160.*1.743847246783397e66-t224+t225+t226+t227-t229.*1.834940238083621e66+t230-t233+t235.*1.955833961275048e66-t238-t239+t241.*1.834940238083621e66-t244+t246.*2.108024684167879e67-t247-t248-t249-t250-t251-t252+t253-t254-t255-t256+t259
	t261 = q1+q2+q4
	t262 = cos(t261)
	t263 = q1+q2-q4+t3
	t264 = cos(t263)
	t265 = q1-q2+q4+t3
	t266 = cos(t265)
	t267 = q1+q2-q4
	t268 = cos(t267)
	t269 = q1-q2+q4
	t270 = cos(t269)
	t271 = -q1+q2+q4
	t272 = cos(t271)
	t273 = q1-q2-q4+t3
	t274 = cos(t273)
	t275 = q1+q2+q4+t3
	t276 = cos(t275)
	t277 = t87.*2.367169127171695e68
	t278 = t164.*2.085484849158977e68
	t279 = t167.*1.222349859852932e66
	t280 = t95.*2.258991771426151e65
	t281 = t99.*1.798910843862435e66
	t282 = t101.*3.838848983032198e67
	t283 = t103.*1.641122594060995e67
	t284 = t272.*1.230830448468477e68
	t285 = t105.*1.798910843862435e66
	t286 = t107.*1.641122594060995e67
	t287 = t274.*9.450565589682945e67
	t288 = t109.*1.777598974740974e67
	t289 = t110.*3.629158332897017e68
	t290 = t112.*1.88536218576757e68
	t291 = t114.*2.305729983768969e67
	t292 = t184.*2.494834337222995e67
	t293 = t118.*1.559651403293754e67
	t294 = t120.*1.559651403293754e67
	t295 = t276.*9.733981529711172e66
	t316 = t174.*8.305405998517524e67
	t317 = t176.*8.305405998517524e67
	t296 = t161+t162+t165+t169+t172+t179+t182+t185+t188+t191-t262.*9.733981529711172e66-t264.*1.0665512025487e67-t266.*1.053915668265494e68+t268.*1.0665512025487e67+t270.*1.367037780294489e68+t277-t278-t279-t280+t281-t282-t283-t284-t285+t286+t287+t288-t289-t290-t291-t292-t293-t294+t295-t316-t317
	t297 = q4+t2
	t298 = cos(t297)
	t299 = q4+t4
	t300 = cos(t299)
	t301 = q4+t3
	t302 = cos(t301)
	t303 = -q4+t2+t3-t4
	t304 = cos(t303)
	t305 = t304.*1.680965282298581e66
	t306 = -q4+t2+t3
	t307 = cos(t306)
	t308 = t307.*9.315304957758298e65
	t309 = t31.*1.798910843862435e66
	t310 = t33.*1.798910843862435e66
	t311 = t19.*1.590481329583487e66
	t312 = -q4+t2
	t313 = cos(t312)
	t314 = -q4+t3
	t315 = cos(t314)
	t318 = cos(q4)
	t319 = t318.*1.704625264978413e69
	t320 = q4-t2+t4
	t321 = cos(t320)
	t322 = q4+t2+t3
	t323 = cos(t322)
	t324 = t323.*9.315304957758298e65
	t325 = q4-t3+t4
	t326 = cos(t325)
	t327 = t326.*1.0665512025487e67
	t328 = q4+t3+t4
	t329 = cos(t328)
	t330 = t329.*1.0665512025487e67
	t331 = t87.*4.16436796838023e82
	t332 = q1+q2-q3-q4
	t333 = cos(t332)
	t334 = t333.*1.67831086311752e81
	t335 = q1-q2+q3-q4
	t336 = cos(t335)
	t337 = t336.*1.462987332903369e82
	t338 = q1-q2-q3+q4
	t339 = cos(t338)
	t340 = t339.*1.658427753522667e82
	t341 = t235.*2.752591593499022e81
	t342 = t99.*3.164677424254534e80
	t343 = t243.*2.9260796499974e82
	t344 = t103.*2.887093399641563e81
	t345 = t272.*2.165299823858942e82
	t346 = t105.*3.164677424254534e80
	t347 = q1+q2+q3+q4
	t348 = cos(t347)
	t349 = t348.*1.531726456607094e81
	t350 = t107.*2.887093399641563e81
	t351 = t274.*1.662561080786584e82
	t352 = t109.*3.127185187600635e81
	t353 = t112.*3.316764233312988e82
	t354 = t114.*4.056283084265199e81
	t355 = q1+q2+q3-q4
	t356 = cos(t355)
	t357 = t356.*1.67831086311752e81
	t358 = q1+q2-q3+q4
	t359 = cos(t358)
	t360 = t359.*1.531726456607094e81
	t361 = q1-q2+q3+q4
	t362 = cos(t361)
	t363 = t362.*1.631976306212276e82
	t364 = -q1+q2+q3+q4
	t365 = cos(t364)
	t366 = t365.*1.487128499210077e82
	t367 = t229.*2.582448803900164e81
	t368 = t232.*2.830742306435468e80
	t369 = t237.*2.797206231116802e81
	t370 = t241.*2.582448803900164e81
	t371 = t246.*2.966780994408209e82
	t372 = t258.*2.830742306435468e80
	t373 = t262.*1.71242014023588e81
	t374 = t264.*1.876296718113235e81
	t375 = t266.*1.854068051125159e82
	t376 = t268.*1.876296718113235e81
	t377 = t270.*2.404918296068613e82
	t378 = t276.*1.71242014023588e81
	t379 = t331+t334+t337+t340+t341+t342-t343-t344-t345-t346+t349+t350+t351+t352-t353-t354-t357-t360-t363-t366-t367-t368-t369+t370+t371+t372-t373-t374-t375+t376+t377+t378
	t380 = t75.*t379.*3.2e11
	t381 = t25.*3.990216591296235e80

	return reshape([0.0,0.0,0.0,0.0,0.0,(t6.*-1.505094268798949e95-t7.*2.138061958130637e94-t9.*1.335878922349698e93+t13.*1.012696775761451e92+t15.*1.012696775761451e92-t17.*2.766549601640899e93-t19.*1.758203636729951e93-t53.*2.365421030580428e93+t55.*3.262055717410504e94+t57.*2.229783408622465e93+t59.*2.229783408622465e93+t61.*2.543398625006958e91+t62.*1.40737488355328e25+t63.*1.40737488355328e25+t64.*1.40737488355328e25-t66.*2.543398625006958e91-t67.*4.084599781729617e94+t69.*2.365421030580428e93-t71.*1.756011369675405e93-t73.*1.756011369675405e93+t138.*1.40737488355328e25+2.076488143314067e95)./(t6.*-3.89576798740676e93-t7.*4.982848297232157e92-t13.*3.170086488882676e89-t15.*3.170086488882676e89-t17.*8.655599649991668e91-t26.*1.025217124508562e92-t27.*5.896511208140987e91+t28+t29-t31.*3.170086488882676e89-t33.*3.170086488882676e89+t34+t37+t40+t43+t46+t49+t50+t51-cos(q2.*-2.0-q4.*2.0+t2+t3).*4.996295833675426e88+4.920249036048467e93),t137,t75.*t223.*-5.62949953421312e25,t75.*t296.*-5.62949953421312e25,t75.*(t62+t63+t64+t122+t123+t124+t129+t130+t134-t139-t140-t382-t383-t384-t385-t386-t387-t388+1.324134645979167e70).*-1.40737488355328e25,0.0,0.0,0.0,0.0,0.0,t137,t75.*(t6.*-1.599827389073318e70-t11.*1.901158272315907e67+t19.*6.36192531833395e66-t26.*1.185068789477605e68+t31.*7.19564337544974e66+t33.*7.19564337544974e66-t48.*1.439128675089948e67+t63+t64+t89.*1.30635276721865e69+t91.*1.30635276721865e69+t93.*2.814005591778681e68-t95.*1.807193417140921e66+t97.*1.648678676226304e67-t101.*3.071079186425758e68-t110.*2.903326666317613e69+t116.*3.366352806800849e67-t118.*1.247721122635003e68-t120.*1.247721122635003e68+t130+t138-t139-t140+1.982648382737728e70).*1.40737488355328e25,t75.*t260.*-4.503599627370496e26,t75.*(t6.*-1.325983974360163e69-t11.*5.28131009027995e66-t26.*2.962671973694012e67-t48.*3.59782168772487e66+t161+t162+t165+t168+t169+t172+t179+t182+t185+t188+t191-t278-t279-t280-t282-t289-t292-t293-t294-t298.*1.73521011604649e67-t300.*2.1331024050974e67-t302.*6.766040747477027e68+t305+t308+t309+t310+t311-t313.*1.73521011604649e67-t315.*6.766040747477027e68-t316-t317+t319-t321.*3.131221120289946e67+t324+t327+t330+1.646284341896401e69).*-5.62949953421312e25,t75.*(t62+t63+t64+t125+t126+t127+t128+t130+t135-t139-t140-t383-t384-t389-t390-t391-t392-t393+1.324134645979167e70).*-1.40737488355328e25,0.0,0.0,0.0,0.0,0.0,t75.*t223.*-5.62949953421312e25,t75.*t260.*-4.503599627370496e26,t75.*(t7.*-9.718123450432042e82-t9.*3.710996046322563e82+t11.*3.192553938622998e80+t13.*2.531741939403627e81+t15.*2.531741939403627e81-t17.*2.969150865079421e82-t19.*4.395509091824878e82-t26.*3.51682660384109e82+t36.*2.025635079020705e81+t39.*2.025635079020705e81-t194.*7.863394511673726e82-t198.*4.475529969786883e82-t200.*4.475529969786883e82+t208.*7.138342018844889e80+t210.*8.368996706703245e83-t214.*6.638827342609521e82+t218.*4.529187690296749e81+t221.*4.529187690296749e81+t381+9.038595841607977e83).*4.0e10,t380,t75.*(t53.*2.902375585010773e60-t55.*4.002547854592237e61-t57.*2.735947994619766e60-t59.*2.735947994619766e60-t61.*3.120754392869666e58+t66.*3.120754392869666e58+t67.*5.011810805674362e61-t69.*2.902375585010773e60+t71.*2.154628905576533e60+t73.*2.154628905576533e60+t142.*3.580201084221302e61+t144.*2.791454010505172e58-t146.*2.596118421008013e60+t148.*1.899359042642789e60+t150.*1.899359042642789e60-t152.*3.531736976202944e61-t154.*2.791454010505172e58+t156.*2.596118421008013e60-t158.*1.92727358274784e60-t160.*1.92727358274784e60).*-4.074974036435136e32,0.0,0.0,0.0,0.0,0.0,t75.*t296.*-5.62949953421312e25,t75.*(t6.*-1.325983974360163e69-t11.*5.28131009027995e66-t26.*2.962671973694012e67-t48.*3.59782168772487e66-t95.*2.258991771426151e65-t101.*3.838848983032198e67-t110.*3.629158332897017e68-t118.*1.559651403293754e67-t120.*1.559651403293754e67+t161+t162-t164.*2.085484849158977e68+t165-t167.*1.222349859852932e66+t168+t169+t172-t174.*8.305405998517524e67-t176.*8.305405998517524e67+t179+t182-t184.*2.494834337222995e67+t185+t188+t191-t298.*1.73521011604649e67-t300.*2.1331024050974e67-t302.*6.766040747477027e68+t305+t308+t309+t310+t311-t313.*1.73521011604649e67-t315.*6.766040747477027e68+t319-t321.*3.131221120289946e67+t324+t327+t330+1.646284341896401e69).*-5.62949953421312e25,t380,t75.*(t6.*-1.257359213783462e85+t9.*1.402621882323949e82-t11.*7.432783173316508e81-t17.*2.612736448493452e83+t19.*1.710517335808237e82-t26.*3.186269261421593e83-t27.*1.779891674226747e83+t31.*2.531741939403627e81+t33.*2.531741939403627e81+t42.*8.899458371133735e82+t45.*8.899458371133735e82-t48.*5.063483878807254e81-t298.*4.884182270022806e82-t300.*6.004149497962351e82-t302.*1.904471161819446e84+t304.*4.731496636864144e81+t307.*2.622025246037675e81-t313.*4.884182270022806e82-t315.*1.904471161819446e84+t318.*4.798093567601945e84-t321.*8.813603919095266e82+t323.*2.622025246037675e81+t326.*3.002074748981175e82+t329.*3.002074748981175e82+t381+1.581779989358623e85).*4.0e10,t75.*(t89.*2.819846300254332e60+t91.*2.819846300254332e60+t93.*6.074211695334503e59-t95.*3.900942991087083e57+t97.*3.558778748073662e58-t101.*6.62912155042827e59-t110.*6.267016968070755e60+t116.*7.266488613040416e58-t118.*2.693286131970666e59-t120.*2.693286131970666e59-t164.*3.601322328061862e60-t167.*2.110816506179317e58+t171.*3.93192726512915e59-t174.*1.434220156437109e60-t176.*1.434220156437109e60+t178.*3.240524874325339e60+t181.*2.312819144102071e58-t184.*4.308208044321151e59+t187.*1.596814857535384e60+t190.*1.596814857535384e60).*3.259979229148109e33],10,4)
end

function Δ(x⁻)
    D_e_num = D_e(x⁻)
    E_2_num = E_2(x⁻)

	tmp = [Diagonal(ones(5)); zeros(2,5)]
    delta_F2 = -(E_2_num*(D_e_num\E_2_num'))\E_2_num*tmp
    delta_bar_dq_e = D_e_num\E_2_num'*delta_F2 + tmp

    R = [0 1 0 0 0; 1 0 0 0 0; 0 0 0 1 0; 0 0 1 0 0; 0 0 0 0 1] # relabelling matrix. Not be confused with a rotation matrix

    delta_dq = [R zeros(5,2)]*delta_bar_dq_e

    x⁺ = [R*x⁻[1:5]; delta_dq*x⁻[6:10]];
end

function D_e(in1)

	#    This function was generated by the Symbolic Math Toolbox version 8.0.
	#    10-Mar-2020 20:42:48

	q1 = in1[1]
	q2 = in1[2]
	q3 = in1[3]
	q4 = in1[4]
	q5 = in1[5]
	t2 = cos(q3);
	t3 = q1-q2+q3-q4;
	t4 = cos(t3);
	t5 = -q1+q2+q4;
	t6 = cos(t5);
	t7 = q1-q2;
	t8 = cos(t7);
	t9 = q1-q2+q3;
	t10 = cos(t9);
	t11 = t2.*1.12675666488;
	t12 = q1+q3+q5;
	t13 = q1+q5;
	t15 = t4.*1.35164722e-2;
	t16 = t10.*8.013743276e-2;
	t18 = t8.*7.1681372485e-2;
	t19 = t6.*1.2090221075e-2;
	t14 = -t15-t16-t18-t19;
	t17 = cos(q4);
	t20 = t17.*2.418044215e-2;
	t21 = q2+q4+q5;
	t22 = q2+q5;
	t23 = t2.*5.6337833244e-1;
	t24 = t23+7.02717990725e-1;
	t25 = -t15-t16;
	t26 = q1+q3;
	t27 = cos(t26);
	t28 = cos(t12);
	t29 = sin(t12);
	t30 = t29.*2.29192575;
	t31 = t17.*1.2090221075e-2;
	t32 = t31+1.0321331925e-2;
	t33 = cos(t21);
	t34 = t33.*4.388465e-2;
	t35 = sin(t21);
	t36 = cos(q1);
	t37 = -t15-t16-t18-t19+t20+7.8872363162e-2;
	t38 = -t15-t16+t23-t27.*5.57372816e-2+7.02717990725e-1;
	t39 = -t15-t19+t31+1.0321331925e-2;
	t40 = cos(t13);
	t41 = cos(t22);
	t42 = t41.*2.6018647e-1;
	t43 = sin(t13);
	t44 = t43.*1.82915043;
	t45 = sin(t22);
	t48 = t28.*2.29192575;
	t49 = t40.*1.82915043;
	t46 = -t48-t49;
	t47 = t34+t42;
	t50 = cos(q5);
	t51 = t50.*1.809652e-1;
	t52 = t30+t44;
	t54 = t35.*4.388465e-2;
	t55 = t45.*2.6018647e-1;
	t53 = -t54-t55;
	t56 = sin(q5);

	return reshape([t11+1.203518592942,t14,t24,-t15-t19,t11-t15-t16-t18-t19-t27.*5.57372816e-2-t36.*4.98559126e-2+1.203518592942,t46,t52,t14,t20+7.8872363162e-2,t25,t32,t37,t47,t53,t24,t25,7.02717990725e-1,-t15,t38,-t48,t30,t4.*(-1.35164722e-2)-t6.*1.2090221075e-2,t32,-t15,1.0321331925e-2,t39,t34,-t54,t4.*(-1.35164722e-2)-t6.*1.2090221075e-2-t8.*7.1681372485e-2-t10.*8.013743276e-2+t11-t27.*5.57372816e-2-t36.*4.98559126e-2+1.203518592942,t37,t38,t39,t4.*(-2.70329444e-2)-t6.*2.418044215e-2-t8.*1.4336274497e-1-t10.*1.6027486552e-1+t11+t20-t27.*1.114745632e-1-t36.*9.97118252e-2+1.313677579864,t34+t42-t48-t49+t51,t30+t44-t54-t55-t56.*1.809652e-1,t46,t47,t28.*(-2.29192575),t34,t28.*(-2.29192575)+t34-t40.*1.82915043+t42+t51,7.5838,0.0,t52,t53,t30,t35.*(-4.388465e-2),t30-t35.*4.388465e-2+t44-t45.*2.6018647e-1-t56.*1.809652e-1,0.0,7.5838],7,7)
end

function E_2(in1)

	#    This function was generated by the Symbolic Math Toolbox version 8.0.
	#    10-Mar-2020 20:42:50

	q1 = in1[1]
	q2 = in1[2]
	q3 = in1[3]
	q4 = in1[4]
	q5 = in1[5]
	t2 = q1+q3+q5;
	t3 = cos(t2);
	t4 = q2+q4+q5;
	t5 = cos(t4);
	t6 = t5.*(7.7e1./2.5e2);
	t7 = q1+q5;
	t8 = cos(t7);
	t9 = q2+q5;
	t10 = cos(t9);
	t11 = t10.*2.755e-1;
	t12 = sin(t2);
	t13 = t12.*(7.7e1./2.5e2);
	t14 = sin(t4);
	t15 = sin(t7);
	t16 = t15.*2.755e-1;
	t17 = sin(t9);

	return reshape([t3.*(-7.7e1./2.5e2)-t8.*2.755e-1,t13+t16,t6+t11,t14.*(-7.7e1./2.5e2)-t17.*2.755e-1,t3.*(-7.7e1./2.5e2),t13,t6,t14.*(-7.7e1./2.5e2),t3.*(-7.7e1./2.5e2)+t6-t8.*2.755e-1+t11,t13-t14.*(7.7e1./2.5e2)+t16-t17.*2.755e-1,1.0,0.0,0.0,1.0],2,7)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
include("../base.jl")
include("fitnessFunctions.jl")

setIndividualType(IntegerGenotype{Int64})
setRandomSeed(9)
setStopCondition(maxIterations = 1000)
setEvaluator(llegar50, -1)
setGenerator(randomIntegerGenerator, 10, 0, 10, popSize = 100)
setSelector(tournamentSelector, 3, samplingWithRep = true, nSelected = 1.0)
setCrossoverOperator(uniformCross, nChildren = 2)
setMutationOperator(uniformMutation, 0, 10, 1)
setReplacementOperator(replaceAllPopulation)
setExperimentSummary(displayFitness=false, printDuringExperiment=true, outputFile="hola.txt")

@time runGenJ(verbose=false)

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using TAMode
using Plots
using OrdinaryDiffEq
using DiffEqDevTools
using IterativeSolvers

params = TAMode.param(ones(15) * 0.5)
solInit = TAMode.getAutocrine(params)
params.gasCur = 100.0
prob = ODEProblem(TAMode.TAMreact, solInit, (0.0, 1000.0), params)

setups = [
    Dict(:alg => Rodas4()),
    Dict(:alg => Rodas42()),
    Dict(:alg => Rodas4P()),
    Dict(:alg => Ros4LStab()),
    Dict(:alg => Rodas5()),
    Dict(:alg => Rosenbrock23()), # slow
    Dict(:alg => TRBDF2(autodiff = false, linsolve = LinSolveGMRES())), # slow
    Dict(:alg => ABDF2(autodiff = false, linsolve = LinSolveGMRES())), # slow
    # Dict(:alg => Exprb43()), # slow
    # Dict(:alg => Exprb32()), # slow
]

test_sol = solve(prob, AutoTsit5(Rodas5()), reltol = 1.0e-16, abstol = 1.0e-16)

tols = 1.0 ./ 10.0 .^ LinRange(9, 12, 4)

wp = WorkPrecisionSet(prob, tols, tols, setups; save_everystep = false, appxsol = test_sol, maxiters = Int(1e6))

plot(wp)


text:
using Graphs, Random
m = 3
backbone = barabasi_albert(1000,3)

u = backbone.fadjlist .|> length |> argmax
dm = length(backbone.fadjlist[u]) - m

cutlink = shuffle(backbone.fadjlist[u])[1:dm]
for v ∈ cutlink
    rem_edge!(backbone, u, v)
end

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
function compute_derivatives(obj::OrdinaryDifferentialEquation, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "computeDerivatives", Vector{jdouble}, (jdouble, Vector{jdouble}), arg0, arg1)
end

function compute_main_state_jacobian(obj::ODEJacobiansProvider, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return jcall(obj, "computeMainStateJacobian", Vector{Vector{jdouble}}, (jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end

function compute_parameter_jacobian(obj::ODEJacobiansProvider, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::JString)
    return jcall(obj, "computeParameterJacobian", Vector{jdouble}, (jdouble, Vector{jdouble}, Vector{jdouble}, JString), arg0, arg1, arg2, arg3)
end

function get_dimension(obj::OrdinaryDifferentialEquation)
    return jcall(obj, "getDimension", jint, ())
end

function get_parameters_names(obj::ODEJacobiansProvider)
    return jcall(obj, "getParametersNames", List, ())
end

function is_supported(obj::ODEJacobiansProvider, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end


text:
using Knet
using Knet: Net, params, axpy!
isdefined(:MNIST) || include("mnist.jl")

@knet function model01(x; winit=Gaussian(0,.1), hidden=32, f1=:relu, o...)
    h    = wf(x; out=hidden, f=f1, winit=winit)
    return wf(h; out=10, f=:soft, winit=winit)
end

function train01(; epochs=1000, lr=0.5, adam=false, batchsize=100, seed=0, target=nothing, o...)
    seed > 0 && setseed(seed)
    global f = compile(:model01; o...)
    global d = minibatch(MNIST.xtrn, MNIST.ytrn, batchsize)
    setp(f, lr=lr, adam=adam)
    global w0,w1,w2,dw,w2_w0,w2_w1
    w0 = w1 = w2 = nothing
    for epoch=1:epochs
        zloss = sloss = iter = 0
        for (x,ygold) in d
            ypred = forw(f, x)
            back(f, ygold, softloss)

            if target != nothing
                # we start at w0 and want to reach w2
                # currently we are at w1 and about to take the step dw
                # compare dw and w2-w1
                # compare w2-w1 and w2-w0
                w1 = weights(f)
                dw = gradients(f)
                if w0 == nothing
                    w0 = deepcopy(w1)
                    w2 = weights(target)
                    w2_w0 = deepcopy(w2)
                    w2_w1 = deepcopy(w2)
                    for i=1:length(w2); axpy!(-1,w0[i],w2_w0[i]); end
                end
                for i=1:length(w2); copy!(w2_w1[i],w2[i]); axpy!(-1,w1[i],w2_w1[i]); end
                @printf("w2-w0=%g,1 w2-w1=%g,%g dw=%g,%g\n",
                        vnorm(w2_w0), vnorm(w2_w1), vcos(w2_w1, w2_w0),
                        lr*vnorm(dw), -vcos(dw, w2_w1))
            end

            update!(f)
            zloss += zeroone(ypred, ygold)
            sloss += softloss(ypred, ygold)
            iter += 1
        end
        println((epoch, sloss/iter, zloss/iter))
        zloss == 0 && break
    end
    return f
end

weights(f::Net)=map(r->r.out, params(f))
gradients(f::Net)=map(r->r.dif, params(f))

function vnorm(w)
    sqnorm = 0
    for u in w; sqnorm += vecnorm(u)^2; end
    sqrt(sqnorm)
end

function vcos(u,v)
    uv = 0
    for i=1:length(u); uv += CUBLAS.dot(u[i],v[i]); end
    uv / (vnorm(u)*vnorm(v))
end

function test01(f, data, loss)
    sumloss = numloss = 0
    for (x,ygold) in data
        ypred = forw(f, x)
        sumloss += loss(ypred, ygold)
        numloss += 1
    end
    sumloss / numloss
end

f1 = train01(seed=1)
f2 = train01(seed=1, target=f1)

# Can get to zloss=0 with hidden=24 using adam=true lr=0.001


text:
include("rootfinding.jl")

f(x) = sin(x)
df(x) = cos(x)

x, it = bisect(f, 0.1, 4.0)

x, it = newton(f, df, 4.0)

g(x) = cos(x)
#anim = cobweb(g, 0, 1, 0.5, 50)
#gif(anim, "anim_fps15.gif", fps = 1)

g1(x) = x - x^3 - 4*x^2 + 10
g2(x) = sqrt(complex(10/x - 4*x))
g3(x) = sqrt(10 - x^3)/2
g4(x) = sqrt(10/(4+x))
g5(x) = x - (x^3 + 4*x^2 - 10)/(3*x^2 + 8*x)
#=
x1, it1 = fixedpoint(g1, 1.5)
x2, it2 = fixedpoint(g2, 1.5)
x3, it3 = fixedpoint(g3, 1.5)
x4, it4 = fixedpoint(g4, 1.5)
x5, it5 = fixedpoint(g5, 1.5)
println((x1,it1))
println((x2,it2))
println((x3,it3))
println((x4,it4))
println((x5,it5))

f(x) = cos(x)-x
df(x) = -sin(x) - 1
#anim = viznewton(f, df, -3, 10, 6, 100)
#gif(anim, "newton_fps1.gif", fps = 1)

x, it = secant(f, 6, 3)
=#

x = brent(x->cos(x)-x, 0.0, 4.0)
println(x)
x1 = bisect(x->cos(x)-x, 0.0, 4.0)
println(x1)


text:
module InstantiateFromURL

# Deps. `using BinaryProvider` needed to instantiate gen_ commands.
using Pkg, BinaryProvider, Suppressor
using HTTP # for activate_github_path
# TOML module
import Pkg.TOML: parsefile
using Markdown

# Code
include("activate.jl")
include("github_project.jl")

# Exports
export activate_github, activate_github_path # default
export packages_to_default_environment, github_project # IJulia-refactor

end


text:
function adversarial_errors(network, X, Y; coeffs = [-0.74, -0.44],
                            δ = 0.02, solver = GLPK.Optimizer, use_gurobi = false,
                            PGD = false, step_factor = 0.1, iterations = 600)
    """ Computes absolute errors with no adversarial perturbations
        Args:
            - network: neural network (NNV model)
            - X: training examples (each column is an example)
            - Y: labels (each column is a label)
            ---------------
            - coeffs: coefficients to get output (linear combination of network outputs)
            - δ: radius for hyperrectangle to check error within
            - solver: solver for mixed integer programs (ignored if PGD true)
            - use_gurobi: whether we are using gurobi (still need to put in solver)
            - PGD: whether or not to use PGD (if not, solve exactly with MIP Verify)
            - step_factor: step factor for PGD (ignored if PGD false)
            - iterations: number of steps for PGD (ignored if PGD false)
    """

    N_samples = size(Y, 2)

    truevals = vec(coeffs' * Y)
    predvals = [coeffs' * compute_output(network, X[:, i]) for i = 1:N_samples]

    minvals = zeros(N_samples)
    maxvals = zeros(N_samples)

    if PGD
        @showprogress "Getting adversial results..." for i = 1:N_samples
            lbs = X[:, i] .- Float32.(δ * ones(size(X, 1)))
            ubs = X[:, i] .+ Float32.(δ * ones(size(X, 1)))
            maxvals[i] = pgd(network, X[:, i], lbs, ubs, coeffs = coeffs,
                        step_size = step_factor, iterations = iterations)
            minvals[i] = -pgd(network, X[:, i], lbs, ubs, coeffs = -coeffs,
                        step_size = step_factor, iterations = iterations)
        end
    else
        @showprogress "Getting adversial results..." for i = 1:N_samples
            input_set = Hyperrectangle(X[:, i], Float32.(δ * ones(size(X, 1))))
            maxvals[i] = optimize_ouput(network, input_set, coeffs = coeffs,
            solver = solver, use_gurobi = use_gurobi)
            minvals[i] = optimize_ouput(network, input_set, coeffs = coeffs,
            solver = solver, use_gurobi = use_gurobi, maximize = false)
        end
    end

    return max.(abs.(truevals .- maxvals), abs.(truevals .- minvals)),
           max.(abs.(predvals .- maxvals), abs.(predvals .- minvals))
end

function optimize_ouput(network, input_set;
    coeffs = [-0.74, -0.44], solver = GLPK.solver, use_gurobi = false, maximize=true,
    obj_threshold=(maximize ? -Inf : Inf))
    """ Function from Chris to get max output
    """
    # Get your bounds
    bounds = NeuralVerification.get_bounds(Ai2z(), network, input_set; before_act=true)

    # Create your model
    if use_gurobi
        model = Model(with_optimizer(solver, OutputFlag = 0))
    else
        model = Model(with_optimizer(solver))
    end
    z = init_vars(model, network, :z, with_input=true)
    δ = init_vars(model, network, :δ, binary=true)
    # get the pre-activation bounds:
    model[:bounds] = bounds
    model[:before_act] = true

    # Add the input constraint
    NeuralVerification.add_set_constraint!(model, input_set, first(z))

    # Encode the network as an MIP
    encode_network!(model, network, BoundedMixedIntegerLP())
    obj = coeffs'*last(z)
    if maximize
        @objective(model, Max, obj)
        if obj_threshold != -Inf
            @constraint(model, obj >= obj_threshold)
        end
    else
        @objective(model, Min, obj)
        if obj_threshold != Inf
            @constraint(model, obj >= obj_threshold)
        end
    end

    # Set lower and upper bounds
    #for the first layer it's special because it has no ẑ
    set_lower_bound.(z[1], low(bounds[1]))
    set_upper_bound.(z[1], high(bounds[1]))
    for i = 2:length(z)-1
        # Set lower and upper bounds for the intermediate layers
        ẑ_i =  _ẑᵢ₊₁(model, i-1)
        z_i = z[i]
        # @constraint(model, ẑ_i .>= low(bounds[i])) These empirically seem to slow it down?
        # @constraint(model, ẑ_i .<= high(bounds[i]))
        z_low = max.(low(bounds[i]), 0.0)
        z_high = max.(high(bounds[i]), 0.0)
        set_lower_bound.(z_i, z_low)
        set_upper_bound.(model[:z][i], z_high)
    end

    # Set lower and upper bounds for the last layer special because
    # it has no ReLU
    set_lower_bound.(z[end], low(bounds[end]))
    set_upper_bound.(z[end], high(bounds[end]))

    optimize!(model)
    if termination_status(model) == OPTIMAL
        return objective_value(model)
    elseif termination_status(model) == INFEASIBLE
        @warn "Infeasible result, did you have an output threshold? If not, then it should never return infeasible"
        return maximize ? -Inf : Inf
    else
        @assert false "Non optimal result"
    end
end

function fgsm(network, x, lbs, ubs; coeffs = [-0.74, -0.44], step_size = 0.1)
    grad_full = NeuralVerification.get_gradient(network, x)
    grad = grad_full' * coeffs
    return clamp.(x + grad * step_size, lbs, ubs)
end

function pgd(network, x, lbs, ubs; coeffs = [-0.74, -0.44], step_size = 0.01,
    iterations = 600)
    cur_x = x
    for i = 1:iterations
        cur_x = fgsm(network, cur_x, lbs, ubs; coeffs = coeffs, step_size = step_size)
    end
    return coeffs' * compute_output(network, cur_x)
end


text:
abstract type AbstractOperationsProblem end

struct GenericOpProblem<:AbstractOperationsProblem end

mutable struct OperationsTemplate
    transmission::Type{<:PM.AbstractPowerModel}
    devices::Dict{Symbol, DeviceModel}
    branches::Dict{Symbol, DeviceModel}
    services::Dict{Symbol, ServiceModel}
end

"""
    OperationsTemplate(::Type{T}) where {T<:PM.AbstractPowerFormulation}

Creates a model reference of the Power Formulation, devices, branches, and services.

# Arguments
- `model::Type{T<:PM.AbstractPowerFormulation}`:
- `devices::Dict{Symbol, DeviceModel}`: device dictionary
- `branches::Dict{Symbol, BranchModel}`: branch dictionary
- `services::Dict{Symbol, ServiceModel}`: service dictionary

# Example
```julia
template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services)
```
"""
function OperationsTemplate(::Type{T}) where {T<:PM.AbstractPowerModel}

    return  OperationsTemplate(T,
                           Dict{Symbol, DeviceModel}(),
                           Dict{Symbol, DeviceModel}(),
                           Dict{Symbol, ServiceModel}())

end

mutable struct OperationsProblem{M<:AbstractOperationsProblem}
    template::OperationsTemplate
    sys::PSY.System
    canonical::Canonical
end

"""
    OperationsProblem(::Type{M},
    template::OperationsTemplate,
    sys::PSY.System;
    optimizer::Union{Nothing, JuMP.OptimizerFactory}=nothing,
    kwargs...) where {M<:AbstractOperationsProblem,
                      T<:PM.AbstractPowerFormulation}

This builds the optimization model and populates the operation model

# Arguments
- `::Type{M} where {M<:AbstractOperationsProblem, T<:PM.AbstractPowerFormulation} = TestOpProblem`:
The abstract operation model type
- `template::OperationsTemplate`: The model reference made up of transmission, devices,
                                          branches, and services.
- `sys::PSY.System`: the system created using Power Systems

# Output
- `op_problem::OperationsProblem`: The operation model contains the model type, model, Power
Systems system, and optimization model.

# Example
```julia
template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationsProblem(TestOpProblem, template, system; optimizer = optimizer)
```

# Accepted Key Words
- `verbose::Bool`: verbose default is true
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model
- `optimizer::union{Nothing, JuMP.OptimizerFactory} = GLPK_optimizer`: The optimizer gets passed
into the optimization model the default is nothing.
- `initial_conditions::DICKDA`: default of Dict{ICKey, Array{InitialCondition}}
- `parameters::Bool`: enable JuMP parameters
- `use_forecast_data::Bool`: if true, forecast collects the time steps in Power Systems,
if false it runs for one time step
- `initial_time::Dates.DateTime`: initial time of forecast
"""
function OperationsProblem(::Type{M},
                        template::OperationsTemplate,
                        sys::PSY.System;
                        optimizer::Union{Nothing, JuMP.OptimizerFactory}=nothing,
                        kwargs...) where {M<:AbstractOperationsProblem}

    op_problem = OperationsProblem{M}(template,
                          sys,
                          Canonical(template.transmission, sys, optimizer; kwargs...))

    build_op_problem!(op_problem; kwargs...)

    return  op_problem

end
"""
    OperationsProblem(op_problem::Type{M},
                    ::Type{T},
                    sys::PSY.System;
                    kwargs...) where {M<:AbstractOperationsProblem,
                                    T<:PM.AbstractPowerFormulation}

This uses the Abstract Power Formulation to build the model reference and
the optimization model and populates the operation model struct.

# Arguments
- `op_problem::Type{M} = where {M<:AbstractOperationsProblem`: Defines the type of the operation model
- `::Type{T} where T<:PM.AbstractPowerFormulation`: The power formulation used for model ref & optimization model
- `sys::PSY.System`: the system created in Power Systems

# Output
- `op_problem::OperationsProblem`: The operation model contains the model type, model, Power
Systems system, and optimization model.

# Example
```julia
template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationsProblem(TestOpProblem, template, system; optimizer = optimizer)
```

# Accepted Key Words
- `verbose::Bool`: verbose default is true
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model
- `optimizer::union{Nothing, JuMP.OptimizerFactory}`: The optimizer gets passed
into the optimization model the default is nothing.
- `initial_conditions::DICKDA`: default of Dict{ICKey, Array{InitialCondition}}
- `parameters::Bool`: enable JuMP parameters
- `use_forecast_data::Bool`: if true, forecast collects the time steps in Power Systems,
if false it runs for one time step
- `initial_time::Dates.DateTime`: initial time of forecast

"""
function OperationsProblem(::Type{M},
                        ::Type{T},
                        sys::PSY.System;
                        kwargs...) where {M<:AbstractOperationsProblem,
                                          T<:PM.AbstractPowerModel}

    optimizer = get(kwargs, :optimizer, nothing)
    return OperationsProblem{M}(OperationsTemplate(T),
                          sys,
                          Canonical(T, sys, optimizer; kwargs...))

end
"""
    OperationsProblem(::Type{T},
                    sys::PSY.System;
                    kwargs...) where {M<:AbstractOperationsProblem,
                                    T<:PM.AbstractPowerFormulation}

This uses the Abstract Power Formulation to build the model reference and
the optimization model and populates the operation model struct.

***Note:*** the abstract operation model is set to the default operation model

# Arguments
- `op_problem::Type{M}`: Defines the type of the operation model
- `::Type{T} where T<:PM.AbstractPowerFormulation`: The power formulation used for model ref & optimization model
- `sys::PSY.System`: the system created in Power Systems

# Output
- `op_problem::OperationsProblem`: The operation model contains the model type, model, Power
Systems system, and optimization model.

# Example
```julia
template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationsProblem(TestOpProblem, template, system; optimizer = optimizer)
```

# Accepted Key Words
- `verbose::Bool`: verbose default is true
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model
- `optimizer::union{Nothing, JuMP.OptimizerFactory}`: The optimizer gets passed
into the optimization model the default is nothing.
- `initial_conditions::DICKDA`: default of Dict{ICKey, Array{InitialCondition}}
- `parameters::Bool`: enable JuMP parameters
- `use_forecast_data::Bool`: if true, forecast collects the time steps in Power Systems,
if false it runs for one time step
- `initial_time::Dates.DateTime`: initial time of forecast

"""
function OperationsProblem(::Type{T},
                        sys::PSY.System;
                        kwargs...) where {T<:PM.AbstractPowerModel}

    return OperationsProblem(GenericOpProblem,
                         T,
                         sys; kwargs...)

end

get_transmission_ref(op_problem::OperationsProblem) = op_problem.template.transmission
get_devices_ref(op_problem::OperationsProblem) = op_problem.template.devices
get_branches_ref(op_problem::OperationsProblem) = op_problem.template.branches
get_services_ref(op_problem::OperationsProblem) = op_problem.template.services
get_system(op_problem::OperationsProblem) = op_problem.sys

function set_transmission_model!(op_problem::OperationsProblem{M},
                               transmission::Type{T}; kwargs...) where {T<:PM.AbstractPowerModel,
                                                                        M<:AbstractOperationsProblem}

    # Reset the canonical
    op_problem.template.transmission = transmission
    op_problem.canonical = Canonical(transmission,
                                        op_problem.sys,
                                        op_problem.canonical.optimizer_factory; kwargs...)

    build_op_problem!(op_problem; kwargs...)

    return
end

function set_devices_template!(op_problem::OperationsProblem{M},
                          devices::Dict{Symbol, DeviceModel}; kwargs...) where M<:AbstractOperationsProblem

    # Reset the canonical
    op_problem.template.devices = devices
    op_problem.canonical = Canonical(op_problem.template.transmission,
                                        op_problem.sys,
                                        op_problem.canonical.optimizer_factory; kwargs...)

    build_op_problem!(op_problem; kwargs...)

    return
end

function set_branches_template!(op_problem::OperationsProblem{M},
                           branches::Dict{Symbol, DeviceModel}; kwargs...) where M<:AbstractOperationsProblem

    # Reset the canonical
    op_problem.template.branches = branches
    op_problem.canonical = Canonical(op_problem.template.transmission,
                                        op_problem.sys,
                                        op_problem.canonical.optimizer_factory; kwargs...)

    build_op_problem!(op_problem; kwargs...)

    return
end

function set_services_template!(op_problem::OperationsProblem{M},
                           services::Dict{Symbol, DeviceModel}; kwargs...) where M<:AbstractOperationsProblem

    # Reset the canonical
    op_problem.template.services = services
    op_problem.canonical = Canonical(op_problem.template.transmission,
                                        op_problem.sys,
                                        op_problem.canonical.optimizer_factory; kwargs...)

    build_op_problem!(op_problem; kwargs...)

    return
end

function set_device_model!(op_problem::OperationsProblem{M},
                           name::Symbol,
                           device::DeviceModel{D, B}; kwargs...) where {D<:PSY.Injection,
                                                                        B<:AbstractDeviceFormulation,
                                                                        M<:AbstractOperationsProblem}

    if haskey(op_problem.template.devices, name)
        op_problem.template.devices[name] = device
        op_problem.canonical = Canonical(op_problem.template.transmission,
                                            op_problem.sys,
                                            op_problem.canonical.optimizer_factory; kwargs...)
        build_op_problem!(op_problem; kwargs...)
    else
        error("Device Model with name $(name) doesn't exist in the model")
    end

    return

end

function set_branch_model!(op_problem::OperationsProblem{M},
                           name::Symbol,
                           branch::DeviceModel{D, B}; kwargs...) where {D<:PSY.Branch,
                                                                        B<:AbstractDeviceFormulation,
                                                                        M<:AbstractOperationsProblem}

    if haskey(op_problem.template.branches, name)
        op_problem.template.branches[name] = branch
        op_problem.canonical = Canonical(op_problem.template.transmission,
                                            op_problem.sys,
                                            op_problem.canonical.optimizer_factory; kwargs...)
        build_op_problem!(op_problem; kwargs...)
    else
        error("Branch Model with name $(name) doesn't exist in the model")
    end

    return

end

function set_services_model!(op_problem::OperationsProblem{M},
                             name::Symbol,
                             service::DeviceModel; kwargs...) where M<:AbstractOperationsProblem

    if haskey(op_problem.template.devices, name)
        op_problem.template.services[name] = service
        op_problem.canonical = Canonical(op_problem.template.transmission,
                                            op_problem.sys,
                                            op_problem.canonical.optimizer_factory; kwargs...)
        build_op_problem!(op_problem; kwargs...)
    else
        error("Branch Model with name $(name) doesn't exist in the model")
    end

    return

end

function construct_device!(op_problem::OperationsProblem,
                           name::Symbol,
                           device_model::DeviceModel;
                           kwargs...)

    if haskey(op_problem.template.devices, name)
        error("Device with model name $(name) already exists in the Opertaion Model")
    end

    devices_ref = get_devices_ref(op_problem)
    devices_ref[name] = device_model

    construct_device!(op_problem.canonical,
                      get_system(op_problem),
                      device_model,
                      get_transmission_ref(op_problem);
                      kwargs...)

    JuMP.@objective(op_problem.canonical.JuMPmodel,
                    MOI.MIN_SENSE,
                    op_problem.canonical.cost_function)

    return

end

function construct_network!(op_problem::OperationsProblem; kwargs...)

    construct_network!(op_problem, op_problem.template.transmission; kwargs...)

    return
end

function construct_network!(op_problem::OperationsProblem,
                            system_formulation::Type{T};
                            kwargs...) where {T<:PM.AbstractPowerModel}

    construct_network!(op_problem.canonical, get_system(op_problem), T; kwargs...)

    return
end

function get_initial_conditions(op_problem::OperationsProblem)
    return op_problem.canonical.initial_conditions
end

function get_initial_conditions(op_problem::OperationsProblem,
                                ic::InitialConditionQuantity,
                                device::PSY.Device)

    canonical = op_problem.canonical
    key = ICKey(ic, device)

    return get_initial_conditions(canonical, key)

end


text:
function nudgeSeaice!(
    ocn :: Ocean;
    τ   :: Float64,
    cfgs...
)
    do_convadjust = cfgs[:do_convadjust]

    Δt = cfgs[:Δt]
    r = Δt / τ
    rr = r / (1.0 + r)

    @loop_hor ocn i j let

        ifrac_clim = ocn.in_flds.IFRACclim[i, j]

        # This provide an option for nudging area. For example, to nudge the seaice in northern hemisphere ONLY,
        # user can simply set IFRAC_clim in the sounthern hemisphere to non-positive values.
        if ifrac_clim >= 0.0

            ifrac = ocn.in_flds.ifrac[i, j]
            T_ML = ocn.T_ML[i, j]
            FLDO = ocn.FLDO[i, j]
            h_ML = ocn.h_ML[i, j]

            # observe that between -1~-2 degC, slope of SST - IFRAC is
            # roughly 100% / 1K. Use this as a diagnostic relation to determine
            # nudged SST.
            ΔT = rr * (ifrac - ifrac_clim)

            T_ML += ΔT
            ocn.T_ML[i, j] = T_ML
            if FLDO > 1
                ocn.Ts[1:FLDO-1, i, j] .= T_ML
            elseif FLDO == -1
                ocn.Ts[1:ocn.Nz[i, j], i, j] .= T_ML
            end

            ocn.seaice_nudge_energy[i, j] = ΔT * ocn.h_ML[i, j] * ρc_sw   / Δt   # + => warming

            OC_updateB!(ocn, i, j)

        else
            ocn.seaice_nudge_energy[i, j] = 0.0
        end

    end

    if do_convadjust
        @loop_hor ocn i j let
            if ocn.in_flds.IFRACclim[i, j] >= 0.0
                OC_doConvectiveAdjustment!(ocn, i, j;)
            end
        end
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using Pkg; Pkg.instantiate()
using BinaryBuilder

name = "HOOMDBlue"
version = v"3.0.0"

# Collection of sources required to complete build
sources = [
    GitSource(
        "https://github.com/glotzerlab/hoomd-blue.git",
        "30b00532c427626b91ea91cd9403a9ed9a4abd7b"
    ),
    GitSource(
        "https://github.com/USCiLab/cereal.git",
        "34eb6f6bd6783018354c7043d5d6aa2eec4e4dbe"
    ),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cereal
mkdir build && cd build

cmake .. \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN} \
    -DJUST_INSTALL_CEREAL=TRUE
cmake --build . --target install -j ${nproc}

cd ..
rm -rf build

cd $WORKSPACE/srcdir/hoomd-blue
git submodule update --init
mkdir build && cd build

cmake .. \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DPYTHON_EXECUTABLE=$prefix/bin/python3 \
    -DENABLE_MPI=ON \
    -DENABLE_GPU=ON \
    -DCMAKE_CUDA_COMPILER=$prefix/cuda/bin/nvcc
cmake --build . --target install -j ${nproc}

cd ..
install_license LICENSE
rm -rf build

LIBNVBLAS_PATH=$(find ${WORKSPACE}/destdir/ -name libnvblas.so)
cd $(dirname ${LIBNVBLAS_PATH})
cat << EOF > nvblas.conf
NVBLAS_CPU_BLAS_LIB $WORKSPACE/lib/libopenblas64_.so
EOF

cd $WORKSPACE/destdir
mv hoomd/include/* include/
rm -rf hoomd/include
mv hoomd lib/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("_md", :libmd, ["lib/hoomd/md"]),
    LibraryProduct("_dem", :libdem, ["lib/hoomd/dem"]),
    LibraryProduct("_mpcd", :libmpcd, ["lib/hoomd/mpcd"]),
    LibraryProduct("_example_plugin", :libplugin, ["lib/hoomd/example_plugin"]),
    LibraryProduct("_metal", :libmetal, ["lib/hoomd/metal"]),
    LibraryProduct("_hpmc", :libhpmc, ["lib/hoomd/hpmc"]),
    LibraryProduct("_hoomd", :libhoomd, ["lib/hoomd"]),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CUDA_full_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("Eigen_jll"),
    Dependency("OpenBLAS_jll"),
    Dependency("OpenMPI_jll"),
    Dependency("pybind11_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: pabloferz/HOOMDBlueBuilder

max_issues_repo_name: pabloferz/HOOMDBlueBuilder

max_forks_repo_name: pabloferz/HOOMDBlueBuilder


text:
module MultiConvex

using Compat, Convex

# package code goes here
include("variablesin.jl")
include("multivexity.jl")
include("partition.jl")
include("solution.jl")

end # module


text:
using Plots, LaTeXStrings
import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum

dir = joinpath(@__DIR__, "src", "assets")

#--------------------------------------------------------------------------------
# ccm89

w = range(0.3, 10.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = ccm89_invum.(w, rv, Ref(ccm89_ca), Ref(ccm89_cb))
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "ccm89_plot.svg"))

#--------------------------------------------------------------------------------
# od94

w = range(0.3, 10.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = ccm89_invum.(w, rv, Ref(od94_ca), Ref(od94_cb))
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "od94_plot.svg"))

#--------------------------------------------------------------------------------
# cal00

w = range(0.46, 8.3, length=1000)
plot()
for rv in [2.0, 3.0, 4.05, 5.0, 6.0]
  m = cal00_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "cal00_plot.svg"))

#--------------------------------------------------------------------------------
# sfd98

l = range(-pi, pi, length=400)
b = range(-pi/64, pi/64, length=300)
dustmap = SFD98Map()
m = [dustmap(li, bj) for li in l, bj in b]
heatmap(l, b, m, label="", transpose=true, colorbar_title="E(B-V)")
xlabel!("l (rad)")
ylabel!("b (rad)")
savefig(joinpath(dir, "sfd98_plot.svg"))


text:
using Test
using JuMP
using MultivariateMoments
using PolyJuMP
using DynamicPolynomials

function _zero_polynomial_test(
    optimizer::MOI.AbstractOptimizer,
    config::MOI.Test.Config,
    complex::Bool,
)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)
    @variable(model, β ≤ 1)
    @variable(model, γ)

    @polyvar x y
    z = complex ? 0 + 0im : 0
    cref = @constraint(model, (α - β) * x * y == z)
    p = γ * x * y
    cγ = @constraint(model, p in PolyJuMP.ZeroPoly())

    @objective(model, Max, α + γ)
    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 1.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 1.0 atol=atol rtol=rtol
    @test value(β) ≈ 1.0 atol=atol rtol=rtol
    @test value(γ) ≈ 0.0 atol=atol rtol=rtol
    @test value(UpperBoundRef(β)) ≈ 1.0 atol=atol rtol=rtol
    @test value(cref) isa MultivariatePolynomials.AbstractPolynomial{Float64}
    @test value(cref) ≈ 0.0 * x * y atol=atol rtol=rtol
    @test value(cγ) ≈ 0.0 * x * y atol=atol rtol=rtol

    @test dual_status(model) == MOI.FEASIBLE_POINT
    @test dual(UpperBoundRef(β)) ≈ -1.0 atol=atol rtol=rtol

    for μ in [dual(cref), moments(cref), dual(cγ), moments(cγ)]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 1
        @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == x*y
    end

    F = MOI.VectorAffineFunction{Float64}
    S = PolyJuMP.ZeroPolynomialSet{FullSpace,MB.MonomialBasis,Monomial{true},
                                   MonomialVector{true}}
    @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
        (MOI.VariableIndex, MOI.LessThan{Float64}), (F, S), (MOI.VectorOfVariables, S)])
    @testset "Delete" begin
        test_delete_bridge(model, cref, 3, ((F, MOI.Zeros, 0),))
    end
end

function real_zero_polynomial_test(
    optimizer::MOI.AbstractOptimizer,
    config::MOI.Test.Config,
)
    return _zero_polynomial_test(optimizer, config, false)
end
function complex_zero_polynomial_test(
    optimizer::MOI.AbstractOptimizer,
    config::MOI.Test.Config,
)
    return _zero_polynomial_test(optimizer, config, false)
end

linear_tests["zero_polynomial"] = real_zero_polynomial_test
linear_tests["zero_polynomial"] = complex_zero_polynomial_test