text:
"""
    run_clust(data::ClustData;norm_op::String="zscore",norm_scope::String="full",method::String="kmeans",representation::String="centroid",n_clust::Int=5,n_init::Int=100,iterations::Int=300,save::String="",attribute_weights::Dict{String,Float64}=Dict{String,Float64}(),get_all_clust_results::Bool=false,kwargs...)
norm_op: "zscore", "01"(not implemented yet)
norm_scope: "full","sequence","hourly"
method: "kmeans","kmedoids","kmedoids_exact","hierarchical"
representation: "centroid","medoid"
"""
function run_clust(data::ClustData;
      norm_op::String="zscore",
      norm_scope::String="full",
      method::String="kmeans",
      representation::String="centroid",
      n_clust::Int=5,
      n_seg::Int=data.T,
      n_init::Int=100,
      iterations::Int=300,
      attribute_weights::Dict{String,Float64}=Dict{String,Float64}(),
      save::String="",#QUESTION dead?
      get_all_clust_results::Bool=false,
      kwargs...
    )

    # When adding new methods: add combination of clust+rep to sup_kw_args
    check_kw_args(norm_op,norm_scope,method,representation)

    # normalize
    # TODO: implement 0-1 normalization and add as a choice to runclust
    data_norm = z_normalize(data;scope=norm_scope)
    if !isempty(attribute_weights)
      data_norm = attribute_weighting(data_norm,attribute_weights)
    end
    data_norm_merged = ClustDataMerged(data_norm)

    #clustering
    b_merged, cost, cost_best, iter =run_clust(data_norm_merged, data; method=method, representation=representation, n_clust=n_clust, n_init=n_init, iterations=iterations, orig_k_ids=deepcopy(data.k_ids), kwargs...)

     if n_seg!=b_merged.T &&  n_seg!=0
       b_merged=intraperiod_segmentation(b_merged;n_seg=n_seg,norm_scope=norm_scope,iterations=iterations)
     else
       n_seg=b_merged.T
     end

    # transfer into ClustData format
    best_results = ClustData(b_merged)
    clust_config = set_clust_config(;norm_op=norm_op, norm_scope=norm_scope, method=method, representation=representation, n_clust=n_clust, n_seg=n_seg, n_init=n_init, iterations=iterations, attribute_weights=attribute_weights)
    # save all locally converged solutions and the best into a struct

    if get_all_clust_results
      clust_result = ClustResultAll(best_results,b_merged.k_ids,cost_best,data_norm_merged.data_type,clust_config,b_merged.centers,b_merged.weights,b_merged.k_ids,cost,iter)
    else
      clust_result =  ClustResultBest(best_results,b_merged.k_ids,cost_best,data_norm_merged.data_type,clust_config)
    end
    #TODO save in save file
    return clust_result
end

"""
function run_clust(data_norm_merged::ClustDataMerged;
                  method::String="kmeans",
                  representation::String="centroid",
                  n_clust::Int=5,
                  n_init::Int=100,
                  iterations::Int=300,
                  orig_k_ids::Array{Int64,1}=Array{Int64,1}(),
                  kwargs...)

method: "kmeans","kmedoids","kmedoids_exact","hierarchical"
representation: "centroid","medoid"
"""
function run_clust(data_norm_merged::ClustDataMerged,
                  data::ClustData;
                  method::String="kmeans",
                  representation::String="centroid",
                  n_clust::Int=5,
                  n_init::Int=100,
                  iterations::Int=300,
                  orig_k_ids::Array{Int64,1}=Array{Int64,1}(),
                  kwargs...)
    # initialize data arrays
    centers = Array{Array{Float64},1}(undef,n_init)
    clustids = Array{Array{Int,1},1}(undef,n_init)
    weights = Array{Array{Float64},1}(undef,n_init)
    cost = Array{Float64,1}(undef,n_init)
    iter = Array{Int,1}(undef,n_init)

    # clustering
    for i = 1:n_init
       # TODO: implement shape based clustering methods
       # function call to the respective function (method + representation)
       fun_name = Symbol("run_clust_"*method*"_"*representation)
       centers[i],weights[i],clustids[i],cost[i],iter[i] =
       @eval $fun_name($data_norm_merged,$n_clust,$iterations;$kwargs...)

       # recalculate centers if medoids is used. Recalculate because medoid is not integrally preserving
      if representation=="medoid"
        centers[i] = resize_medoids(data,centers[i],weights[i])
      end
    end
    # find best
    # TODO: write as function
    cost_best,ind_mincost = findmin(cost)  # along dimension 2, only store indice

    k_ids=orig_k_ids
    k_ids[findall(orig_k_ids.!=0)]=clustids[ind_mincost]
    # save in merged format as array
    # NOTE if you need clustered data more precise than 8 digits change the following line accordingly
     n_digits_data_round=8 # Gurobi throws warning when rounding errors on order~1e-13 are passed in. Rounding errors occur in clustering of many zeros (e.g. solar).
     return ClustDataMerged(data_norm_merged.region,data_norm_merged.years,n_clust,data_norm_merged.T,round.(centers[ind_mincost]; digits=n_digits_data_round),data_norm_merged.data_type,weights[ind_mincost],k_ids), cost, cost_best, iter
 end

"""
    run_clust(data::ClustData,n_clust_ar::Array{Int,1};norm_op::String="zscore",norm_scope::String="full",method::String="kmeans",representation::String="centroid",n_init::Int=100,iterations::Int=300,save::String="",kwargs...)
This function is a wrapper function around run_clust(). It runs multiple number of clusters k and returns an array of results.
norm_op: "zscore", "01"(not implemented yet)
norm_scope: "full","sequence","hourly"
method: "kmeans","kmedoids","kmedoids_exact","hierarchical"
representation: "centroid","medoid"
"""
function run_clust(
      data::ClustData,
      n_clust_ar::Array{Int,1};
      norm_op::String="zscore",
      norm_scope::String="full",
      method::String="kmeans",
      representation::String="centroid",
      n_init::Int=100,
      iterations::Int=300,
      save::String="",
      kwargs...
    )
    results_ar = Array{ClustResult,1}(undef,length(n_clust_ar))
    for i=1:length(n_clust_ar)
      results_ar[i] = run_clust(data;norm_op=norm_op,norm_scope=norm_scope,method=method,representation=representation,n_init=n_init,n_clust=n_clust_ar[i],iterations=iterations,save=save,kwargs...)
    end
    return results_ar
end

# supported keyword arguments
sup_kw_args =Dict{String,Array{String}}()
sup_kw_args["region"]=["GER","CA"]
sup_kw_args["opt_problems"]=["battery","gas_turbine"]
sup_kw_args["norm_op"]=["zscore"]
sup_kw_args["norm_scope"]=["full","hourly","sequence"]
sup_kw_args["method+representation"]=["kmeans+centroid","kmeans+medoid","kmedoids+medoid","kmedoids_exact+medoid","hierarchical+centroid","hierarchical+medoid","predefined+centroid","predefined+medoid"]#["dbaclust+centroid","kshape+centroid"]

"""
    get_sup_kw_args
Returns supported keyword arguments for clustering function run_clust()
"""
function get_sup_kw_args()
    return sup_kw_args
end

"""
    check_kw_args(region,opt_problems,norm_op,norm_scope,method,representation)
checks if the arguments supplied for run_clust are supported
"""
function check_kw_args(
      norm_op::String,
      norm_scope::String,
      method::String,
      representation::String
    )
    check_ok = true
    error_string = "The following keyword arguments / combinations are not currently supported: \n"
    # norm_op
    if !(norm_op in sup_kw_args["norm_op"])
       check_ok=false
       error_string = error_string * "normalization operation $norm_op is not supported \n"
    end
    # norm_scope
    if !(norm_scope in sup_kw_args["norm_scope"])
       check_ok=false
       error_string = error_string * "normalization scope $norm_scope is not supported \n"
    end
    # method +  representation
    if !(method*"+"*representation in sup_kw_args["method+representation"])
       check_ok=false
       error_string = error_string * "the combination of method $method and representation $representation is not supported \n"
    elseif method == "dbaclust"
       @info("dbaclust can be run in parallel using src/clust_algorithms/runfiles/cluster_gen_dbaclust_parallel.jl")
    elseif method =="kshape"
       check_ok=false
       error_string = error_string * "kshape is implemented in python and should be run individually: src/clust_algorithms/runfiles/cluster_gen_kshape.py \n"
    end
    error_string = error_string * "get_sup_kw_args() provides a list of supported keyword arguments."

    if check_ok
       return true
    else
       error(error_string)
    end
end

"""
      run_clust_predefined_centroid(data_norm::ClustDataMerged,n_clust::Int,iterations::Int;k_ids::Array{Int64,1}=Array{Int64,1}())
Provide the assigned k_ids before running any clustering within this run
"""
function run_clust_predefined_centroid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int;
    k_ids::Array{Int64,1}=Array{Int64,1}())
    centers,weights,clustids,cost,iter =[],[],[],0,0
    # if only one cluster
    if n_clust == length(unique(k_ids))
        centers_norm=zeros(size(data_norm.data,1),n_clust)
        for i in 1:n_clust
          index=findall(k_ids.==i)
          centers_norm[:,i] = mean(data_norm.data[:,index],dims=2) # should be 0 due to normalization
        end
        clustids = k_ids
        centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids) # need to provide idx in case that sequence-based normalization is used
        cost = sum(pairwise(SqEuclidean(),centers_norm,data_norm.data; dims=2)) #same as sum((seq_norm-repmat(mean(seq_norm,2),1,size(seq,2))).^2)
        iter = 1
    # kmeans() in Clustering.jl is implemented for k>=2
    else
        throw(@error "mismatch of intendet clusters $n_clust and k_ids provided with $(length(unique(k_ids)))")
    end
    weights = calc_weights(clustids,n_clust)
    return centers,weights,clustids,cost,iter
end

"""
      run_clust_predefined_medoid(data_norm::ClustDataMerged,n_clust::Int,iterations::Int;k_ids::Array{Int64,1}=Array{Int64,1}())
Provide the assigned k_ids before running any clustering within this run
"""
function run_clust_predefined_medoid(data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int;
    k_ids::Array{Int64,1}=Array{Int64,1}()
    )
    return run_clust_predefined_centroid(data_norm, n_clust, iterations; k_ids=k_ids)
  end

"""
    run_clust_kmeans_centroid(data_norm::ClustDataMerged,n_clust::Int,iterations::Int)
"""
function run_clust_kmeans_centroid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int
    )
    centers,weights,clustids,cost,iter =[],[],[],0,0
    # if only one cluster
    if n_clust ==1
        centers_norm = mean(data_norm.data,dims=2) # should be 0 due to normalization
        clustids = ones(Int,size(data_norm.data,2))
        centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids) # need to provide idx in case that sequence-based normalization is used
        cost = sum(pairwise(SqEuclidean(),centers_norm,data_norm.data; dims=2)) #same as sum((seq_norm-repmat(mean(seq_norm,2),1,size(seq,2))).^2)
        iter = 1
    # kmeans() in Clustering.jl is implemented for k>=2
    elseif n_clust==data_norm.K
        clustids = collect(1:data_norm.K)
        centers = undo_z_normalize(data_norm.data,data_norm.mean,data_norm.sdv;idx=clustids) # need to provide idx in case that sequence-based normalization is used
        cost = 0.0
        iter = 1
    else
        results = kmeans(data_norm.data,n_clust;maxiter=iterations)
        # save clustering results
        clustids = results.assignments
        centers_norm = results.centers
        centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids)
        cost = results.totalcost
        iter = results.iterations
    end

    weights = calc_weights(clustids,n_clust)

    return centers,weights,clustids,cost,iter

end

"""
    run_clust_kmeans_medoid(
      data_norm::ClustDataMerged,
      n_clust::Int,
      iterations::Int
    )
"""
function run_clust_kmeans_medoid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int
    )
    centers,weights,clustids,cost,iter =[],[],[],0,0
    # if only one cluster
    if n_clust ==1
        clustids = ones(Int,size(data_norm.data,2))
        centers_norm = calc_medoids(data_norm.data,clustids)
        centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids) # need to provide idx in case that sequence-based normalization is used
        cost = sum(pairwise(SqEuclidean(),centers_norm,data_norm.data; dims=2)) #same as sum((seq_norm-repmat(mean(seq_norm,2),1,size(seq,2))).^2)
        iter = 1
    # kmeans() in Clustering.jl is implemented for k>=2
    elseif n_clust==data_norm.K
        clustids = collect(1:data_norm.K)
        centers = undo_z_normalize(data_norm,data_norm.mean,data_norm.sdv;idx=clustids) # need to provide idx in case that sequence-based normalization is used
        cost = 0.0
        iter = 1
    else
        results = kmeans(data_norm.data,n_clust;maxiter=iterations)

        # save clustering results
        clustids = results.assignments
        centers_norm = calc_medoids(data_norm.data,clustids)
        centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids)
        cost = calc_SSE(data_norm.data,centers_norm,clustids)
        iter = results.iterations
    end

    weights = calc_weights(clustids,n_clust)

    return centers,weights,clustids,cost,iter

end

"""
    run_clust_kmedoids_medoid(
      data_norm::ClustDataMerged,
      n_clust::Int,
      iterations::Int
    )
"""
function run_clust_kmedoids_medoid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int
    )

    # TODO: optional in future: pass distance metric as kwargs
    dist = SqEuclidean()
    d_mat=pairwise(dist,data_norm.data; dims=2)
    results = kmedoids(d_mat,n_clust;tol=1e-6,maxiter=iterations)
    clustids = results.assignments
    centers_norm = data_norm.data[:,results.medoids]
    centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids)
    cost = results.totalcost
    iter = results.iterations

    weights = calc_weights(clustids,n_clust)

    return centers,weights,clustids,cost,iter
end

"""
    run_clust_kmedoids_exact_medoid(
      data_norm::ClustDataMerged,
      n_clust::Int,
      iterations::Int;
      gurobi_env=0
    )
"""
function run_clust_kmedoids_exact_medoid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int;
    gurobi_env=0
    )

    (typeof(gurobi_env)==Int) && @error("Please provide a gurobi_env (Gurobi Environment). See test file for example")

    # TODO: optional in future: pass distance metric as kwargs
    dist = SqEuclidean()
    results = kmedoids_exact(data_norm.data,n_clust,gurobi_env;_dist=dist)#;distance_type_ar[dist])
    clustids = results.assignments
    centers_norm = results.medoids
    centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids)
    cost = results.totalcost
    iter = 1

    weights = calc_weights(clustids,n_clust)

    return centers,weights,clustids,cost,iter
end

"""
    run_clust_hierarchical(
      data_norm::ClustDataMerged,
      n_clust::Int,
      iterations::Int;
      _dist::SemiMetric = SqEuclidean()
    )

Helper function to run run_clust_hierarchical_centroids and run_clust_hierarchical_medoid
"""
function run_clust_hierarchical(
    data::Array{Float64,2},
    n_clust::Int,
    iterations::Int;
    _dist::SemiMetric = SqEuclidean()
    )

    d_mat=pairwise(_dist,data; dims=2)
    r=hclust(d_mat,linkage=:ward_presquared)
    clustids = cutree(r,k=n_clust)
    weights = calc_weights(clustids,n_clust)

    return [],weights,clustids,[],1
end

"""
    run_clust_hierarchical_centroid(
      data_norm::ClustDataMerged,
      n_clust::Int,
      iterations::Int;
      _dist::SemiMetric = SqEuclidean()
    )
"""
function run_clust_hierarchical_centroid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int;
    _dist::SemiMetric = SqEuclidean()
    )
    x,weights,clustids,x,iter= run_clust_hierarchical(data_norm.data,n_clust,iterations;_dist=_dist)
    centers_norm = calc_centroids(data_norm.data,clustids)
    cost = calc_SSE(data_norm.data,centers_norm,clustids)
    centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids)

    return centers,weights,clustids,cost,iter
end

"""
    run_clust_hierarchical_medoid(
      data_norm::ClustDataMerged,
      n_clust::Int,
      iterations::Int;
      _dist::SemiMetric = SqEuclidean()
    )
"""
function run_clust_hierarchical_medoid(
    data_norm::ClustDataMerged,
    n_clust::Int,
    iterations::Int;
    _dist::SemiMetric = SqEuclidean()
    )
    ~,weights,clustids,~,iter= run_clust_hierarchical(data_norm,n_clust,iterations;_dist=_dist)
    centers_norm = calc_medoids(data_norm.data,clustids)
    cost = calc_SSE(data_norm.data,centers_norm,clustids)
    centers = undo_z_normalize(centers_norm,data_norm.mean,data_norm.sdv;idx=clustids)

    return centers,weights,clustids,cost,iter
end


text:
module AnyMODResults

using CSV
using DataFrames
using Glob
using PlotlyJS
using Pkg

include("postprocessing.jl")
include("AnymodResult.jl")
include("Mask.jl")
include("plots.jl")

end # module

ext: jl


text:
# using Bonsai, HTTP

# function timer_middleware(stream, next)
#     x = now()
#     # the next middleware or handler in the stack
#     next(stream)
#     elapsed = now() - x
#     @info "$(stream.message.target) took $elapsed"
# end

# function error_middleware(stream, next)
#     try
#         next(stream)
#     catch e
#         @error e
#         code = InternalServerError()
#         Bonsai.write(stream, HTTP.statustext(code), code)
#     end
# end

# function error_handler(stream)
#     error("Oh no")
# end

# function index_handler(stream)
#     Bonsai.write(stream, "ok")
# end

# router = Router()
# # make sure error handler is the first middleware
# all!(router, "*", error_middleware)
# all!(router, "*", timer_middleware)
# get!(router, "*", file_handler)
# get!(router, "/error", error_handler)
# wait(start(router))

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
include("hopf.jl")
include("planar_contact.jl")


text:
## --- Read ESRI Arc/Info ASCII grid files

    function importAAIGrid(fname, T=Float64; undefval=NaN)
        # Open the file
        fid = open(fname)

        metadata = Dict{String,Number}()
        metadata["ncols"] = parse(Int64, match(r"  *(.*?)$", readline(fid))[1])
        metadata["nrows"] = parse(Int64, match(r"  *(.*?)$", readline(fid))[1])
        metadata["xll_corner"] = parse(Float64, match(r"  *(.*?)$", readline(fid))[1])
        metadata["yll_corner"] = parse(Float64, match(r"  *(.*?)$", readline(fid))[1])
        metadata["cellsize"] = parse(Float64, match(r"  *(.*?)$", readline(fid))[1])
        metadata["nodata"] = parse(Float64, match(r"  *(.*?)$", readline(fid))[1])

        nrows = metadata["nrows"]
        ncols = metadata["ncols"]

        data = Array{T}(undef,ncols,nrows)
        for i = 1:nrows
            l = readline(fid)
            delim_string_parse!(data, l, ' ', T, offset=(i-1)*ncols, undefval=undefval)
        end

        # Close the file
        close(fid)

        return (data', metadata)
    end
    export importAAIGrid

## --- Calculate slope from a DEM

    function maxslope(matrix, x_lon_cntr, y_lat_cntr, cellsize, T=UInt16; minmatval=-12000, km_per_lat=111.1)
        # Returns slope in units/kilometer given a latitude-longitude grid of z-values

        # Allocate output array
        slope = Array{T}(undef,size(matrix))

        # Fill in the center first
        distNS = 2 * cellsize * km_per_lat
        for i = 2:(size(matrix,1)-1)
            # Distance between grid cell centers
            km_per_lon = cos(y_lat_cntr[i] * pi/180) * km_per_lat
            distEW = 2 * cellsize * km_per_lon
            distDiag = sqrt(distNS^2 + distEW^2)

            for j = 2:(size(matrix,2)-1)
                # Gradients, in matrix units per km
                if matrix[i,j] < minmatval
                    slope[i,j] = 0
                else
                    if (matrix[i+1,j] < minmatval) || (matrix[i-1,j] < minmatval)
                        NS = 0
                    else
                        NS = abs(matrix[i+1,j] - matrix[i-1,j]) / distNS
                    end
                    if (matrix[i,j+1] < minmatval) || (matrix[i,j-1]<minmatval)
                        EW = 0
                    else
                        EW = abs(matrix[i,j+1] - matrix[i,j-1]) / distEW
                    end
                    if (matrix[i+1,j-1] < minmatval) || (matrix[i-1,j+1] < minmatval)
                        NESW = 0
                    else
                        NESW = abs(matrix[i+1,j-1] - matrix[i-1,j+1]) / distDiag
                    end
                    if (matrix[i+1,j+1] < minmatval) || (matrix[i-1,j-1] < minmatval)
                        NWSE = 0
                    else
                        NWSE = abs(matrix[i+1,j+1] - matrix[i-1,j-1]) / distDiag
                    end

                    # Record the steepest slope
                    slope[i,j] = nearest(T, max(NS,EW,NESW,NWSE))
                end
            end

            # Fill in edges too
            distEW = cellsize * km_per_lon
            distDiag = sqrt((distNS/2)^2 + distEW^2)

            # Left edge
            if (matrix[i+1,1] < minmatval) || (matrix[i-1,1] < minmatval)
                NS = 0
            else
                NS = abs(matrix[i+1,1] - matrix[i-1,1]) / distNS
            end
            if (matrix[i,2] < minmatval) || (matrix[i,1] < minmatval)
                EW = 0
            else
                EW = abs(matrix[i,2] - matrix[i,1]) / distEW
            end
            if (matrix[i+1,1] < minmatval) || (matrix[i-1,2] < minmatval)
                NESW = 0
            else
                NESW = abs(matrix[i+1,1] - matrix[i-1,2]) / distDiag
            end
            if (matrix[i+1,2] < minmatval) || (matrix[i-1,1] < minmatval)
                NWSE = 0
            else
                NWSE = abs(matrix[i+1,2] - matrix[i-1,1]) / distDiag
            end
            slope[i,1] = nearest(T, max(NS,EW,NESW,NWSE))

            # Right edge
            if (matrix[i+1,end] < minmatval) || (matrix[i-1,end] < minmatval)
                NS = 0
            else
                NS = abs(matrix[i+1,end] - matrix[i-1,end]) / distNS
            end
            if matrix[i,end]<minmatval || matrix[i,end-1]<minmatval
                EW = 0
            else
                EW = abs(matrix[i,end] - matrix[i,end-1]) / distEW
            end
            if (matrix[i+1,end-1] < minmatval) || (matrix[i-1,end] < minmatval)
                NEWS = 0
            else
                NESW = abs(matrix[i+1,end-1] - matrix[i-1,end]) / distDiag
            end
            if (matrix[i+1,end] < minmatval) || (matrix[i-1,end-1] < minmatval)
                NWSE = 0
            else
                NWSE = abs(matrix[i+1,end] - matrix[i-1,end-1]) / distDiag
            end
            slope[i,end] = nearest(T, max(NS,EW,NESW,NWSE))
        end

        # Fill in the top and bottom row
        distNS = cellsize * km_per_lat

        # Top row
        km_per_lon = cos(y_lat_cntr[1]*pi/180) * km_per_lat
        distEW = 2*cellsize*km_per_lon
        distDiag = sqrt(distNS^2+(distEW/2)^2)

        for j = 2:(size(matrix,2)-1)
            # Gradients, in meters per km
            if (matrix[2,j] < minmatval) || (matrix[1,j] < minmatval)
                NS = 0
            else
                NS = abs(matrix[2,j] - matrix[1,j]) / distNS
            end
            if (matrix[1,j+1] < minmatval) || (matrix[1,j-1] < minmatval)
                EW = 0
            else
                EW = abs(matrix[1,j+1] - matrix[1,j-1]) / distEW
            end
            if (matrix[2,j-1] < minmatval) || (matrix[1,j] < minmatval)
                NESW = 0
            else
                NESW = abs(matrix[2,j-1] - matrix[1,j]) / distDiag
            end
            if (matrix[2,j+1] < minmatval) || (matrix[1,j] < minmatval)
                NWSE = 0
            else
                NWSE = abs(matrix[2,j+1] - matrix[1,j]) / distDiag
            end
            slope[1,j] = nearest(T, max(NS,EW,NESW,NWSE))
        end
        slope[1,1] = 0
        slope[1,end] = 0

        # Bottom row
        km_per_lon = cos(y_lat_cntr[end] * pi/180) * km_per_lat
        distEW = 2 * cellsize * km_per_lon
        distDiag = sqrt(distNS^2 + (distEW/2)^2)
        for j = 2:(size(matrix,2)-1)
            # Gradients, in meters per Km
            if (matrix[end-1,j] < minmatval) || (matrix[end,j] < minmatval)
                NS = 0
            else
                NS = abs(matrix[end-1,j] - matrix[end,j]) / distNS
            end
            if (matrix[end,j+1] < minmatval) || (matrix[end,j-1] < minmatval)
                EW = 0
            else
                EW = abs(matrix[end,j+1] - matrix[end,j-1]) / distEW
            end
            if (matrix[end-1,j-1] < minmatval) || (matrix[end,j] < minmatval)
                NESW = 0
            else
                NESW = abs(matrix[end-1,j-1] - matrix[end,j]) / distDiag
            end
            if (matrix[end-1,j+1] < minmatval) || (matrix[end,j] < minmatval)
                NWSE = 0
            else
                NWSE = abs(matrix[end-1,j+1] - matrix[end,j]) / distDiag
            end
            slope[end,j] = nearest(T, max(NS,EW,NESW,NWSE))
        end
        slope[end,1] = 0
        slope[end,end] = 0

        return slope
    end
    export maxslope

    function aveslope(matrix, x_lon_cntr, y_lat_cntr, cellsize, T=UInt16; minmatval=-12000, maxmatval=9000, km_per_lat=111.1)
        # Returns slope in units/kilometer given a latitude-longitude grid of z-values

        # Allocate intermediate and output arrays
        distance = Array{Float64}(undef,8)
        local_slopes = Array{Float64}(undef,8)
        slope = Array{T}(undef,size(matrix))

        # Index offsets to cycle through:
        #         [N,NE,E,SE,S,SW,W,NW]
        ioffset = [-1,-1,0,1,1,1,0,-1]
        joffset = [0,1,1,1,0,-1,-1,-1]
        #
        # i.e. Layout:
        # 8 1 2
        # 7 x 3
        # 6 5 4

        # Distance between grid cell centers
        # N, S
        distance[[1,5]] .= cellsize * km_per_lat

        # Fill in the center first
        for i = 2:(size(matrix,1)-1)
            # Distance between grid cell centers
            km_per_lon = cos(y_lat_cntr[i]*pi/180) * km_per_lat
            distance[[3,7]] .= cellsize*km_per_lon; #E, W
            distance[[2,4,6,8]] .= sqrt(distance[1]^2+distance[3]^2)  # Diagonals

            # Center
            for j = 2:(size(matrix,2)-1)
                # Gradients, in matrix z-units per km
                here = matrix[i,j]
                if (here < minmatval) || (here > maxmatval)
                    slope[i,j] = 0
                else
                    for k = 1:8
                        there = matrix[i+ioffset[k], j+joffset[k]]
                        if (there < minmatval) || (there > maxmatval)
                            local_slopes[k] = 0
                        else
                            local_slopes[k] = abs(there-here) / distance[k]
                        end
                    end
                    # Record the average slope
                    slope[i,j] = nearest(T, mean(local_slopes))
                end
            end

            # Left edge
            here = matrix[i,1]
            if (here < minmatval) || (here > maxmatval)
                slope[i,1] = 0
            else
                for k = 1:5
                    there = matrix[i+ioffset[k], 1+joffset[k]]
                    if (there < minmatval) || (there > maxmatval)
                        local_slopes[k] = 0
                    else
                        local_slopes[k] = abs(there-here) / distance[k]
                    end
                end
                slope[i,1] = nearest(T, mean(local_slopes[1:5]))
            end

            # Right edge
            here = matrix[i,end]
            if (here < minmatval) || (here > maxmatval)
                slope[i,end] = 0
            else
                for k = [5,6,7,8,1]
                    there = matrix[i+ioffset[k], end+joffset[k]]
                    if (there < minmatval) || (there > maxmatval)
                        local_slopes[k] = 0
                    else
                        local_slopes[k] = abs(there-here) / distance[k]
                    end
                end
                slope[i,end] = nearest(T, mean(local_slopes[[5,6,7,8,1]]))
            end
        end

        # Top row
        km_per_lon = cos(y_lat_cntr[1] * pi/180) * km_per_lat
        distance[[3,7]] .= cellsize * km_per_lon #E, W
        distance[[2,4,6,8]] .= sqrt(distance[1]^2 + distance[3]^2)  # Diagonals
        for j = 2:(size(matrix,2)-1)
            # Gradients, in matrix units per km
            here = matrix[1,j]
            if (here < minmatval) || (here > maxmatval)
                slope[1,j] = 0
            else
                for k=3:7
                    there = matrix[1+ioffset[k], j+joffset[k]]
                    if (there < minmatval) || (there > maxmatval)
                        local_slopes[k] = 0
                    else
                        local_slopes[k] = abs(there-here) / distance[k]
                    end
                end
                slope[1,j] = nearest(T, mean(local_slopes[3:7]))
            end
        end
        slope[1,1] = 0
        slope[1,end] = 0

        # Bottom row
        km_per_lon = cos(y_lat_cntr[end] *pi/180) * km_per_lat
        distance[[3,7]] .= cellsize * km_per_lon #E, W
        distance[[2,4,6,8]] .= sqrt(distance[1]^2+distance[3]^2)  # Diagonals
        for j = 2:(size(matrix,2)-1)
            # Gradients, in matrix units per km
            here = matrix[end,j]
            if (here < minmatval) || (here > maxmatval)
                slope[end,j] = 0
            else
                for k = [7,8,1,2,3]
                    there = matrix[end+ioffset[k], j+joffset[k]]
                    if (there < minmatval) || (there > maxmatval)
                        local_slopes[k] = 0
                    else
                        local_slopes[k] = abs(there-here) / distance[k]
                    end
                end
                slope[end,j] = nearest(T, mean(local_slopes[[7,8,1,2,3]]))
            end
        end
        slope[end,1] = 0
        slope[end,end] = 0

        return slope
    end
    export aveslope

## --- Generate random latitude and longitude pairs uniformly distributed across the globe

    CONST_180_PI = 180/pi

    function randlatlon(n::Integer)
        90 .- CONST_180_PI*acos.(2*rand(n) .- 1), rand(n)*360 .- 180
    end
    function randlatlon()
        90 .- CONST_180_PI*acos.(2*rand() .- 1), rand()*360 .- 180
    end
    export randlatlon

## --- End of File


text:
@discord_object struct GuildMember
    user::User
    nick::String
    roles::Vector{Snowflake}
    joined_at::Union{String, DateTime}
    premium_since::Union{String, DateTime}
    deaf::Bool
    mute::Bool
end


text:
using Flux
using Flux:glorot_uniform
using AtomicGraphNets

"""
Build a model of the architecture introduced in the Xie and Grossman 2018 paper: https://arxiv.org/abs/1710.10324

Input to the resulting model is a FeaturedGraph with feature matrix with `input_feature_length` rows and one column for each node in the input graph.

Network has convolution layers, then pooling to some fixed length, followed by Dense layers leading to output.

# Arguments
- `input_feature_length::Integer`: length of feature vector at each node
- `num_conv::Integer`: number of convolutional layers
- `conv_activation::F`: activation function on convolutional layers
- `atom_conv_feature_length::Integer`: length of output of conv layers
- `pool_type::String`: type of pooling after convolution (mean or max)
- `pool_width::Float`: fraction of atom_conv_feature_length that pooling window should span
- `pooled_feature_length::Integer`: feature length to pool down to
- `num_hidden_layers::Integer`: how many Dense layers before output? Note that if this is set to 1 there will be no nonlinearity imposed on these layers
- `hidden_layer_activation::F`: activation function on hidden layers
- `output_layer_activation::F`: activation function on output layer; should generally be identity for regression and something that normalizes appropriately (e.g. softmax) for classification
- `output_length::Integer`: length of output vector
- `initW::F`: function to use to initialize weights in trainable layers
"""
function Xie_model(input_feature_length; num_conv=2, conv_activation=softplus, atom_conv_feature_length=80, pool_type="mean", pool_width=0.1, pooled_feature_length=40, num_hidden_layers=1, hidden_layer_activation=softplus, output_layer_activation=identity, output_length=1, initW=glorot_uniform)
    @assert atom_conv_feature_length >= pooled_feature_length "Feature length after pooling must be <= feature length before pooling!"
    model = Chain(AGNConv(input_feature_length=>atom_conv_feature_length, conv_activation, initW=initW), [AGNConv(atom_conv_feature_length=>atom_conv_feature_length, conv_activation, initW=initW) for i in 1:num_conv-1]..., AGNPool(pool_type, atom_conv_feature_length, pooled_feature_length, pool_width), [Dense(pooled_feature_length, pooled_feature_length, hidden_layer_activation, initW=initW) for i in 1:num_hidden_layers-1]..., Dense(pooled_feature_length, output_length, output_layer_activation, initW=initW))
end


text:
#ConventionalApp: base type, constant & function definitions.
#-------------------------------------------------------------------------------

#==Main types
===============================================================================#
abstract type AbstractProject; end

struct ActiveProject <: AbstractProject; end
const activeproject = ActiveProject()

struct Project <: AbstractProject
	root::String
end
Project(m::Module) = Project(abspath(dirname(dirname(pathof(m)))))
Project(::ActiveProject) = Project(abspath(dirname(Base.load_path_expand("@"))))
#Project(repo::LibGit2.GitRepo) = Project(LibGit2.path(repo))

#==Extensions
===============================================================================#
Base.basename(proj::Project) = basename(proj.root)
Pkg.activate(p::Project) = Pkg.activate(p.root)

#==Display functions
===============================================================================#

"""
`displayenv()`

Display environment info.
"""
function displayenv()
	projpath = dirname(Base.load_path_expand("@"))
	@info("Project: $(basename(projpath))")
	map((path)->println("--> $path"), Base.load_path())
	println()
end
Base.show(::ActiveProject) = displayenv()

#==Top-level interface
===============================================================================#

"""
`@include_startup(proj::Project)`

`include` the `./startup.jl` file found in the `proj` directory. Should be
called from a global scope (not a function).

Implemented as macro to ensure `include` is run from calling scope.
"""
macro include_startup(proj)
	return quote begin #Create temp scope for startupfile variable
		local startupfile = joinpath($(esc(proj)).root, "startup.jl")
		if ispath(startupfile)
			$(__module__).include(startupfile)
		end
	end end
end

"""
`setup_env(proj::Project)`

Define a minmial environment stack and activates specified `Project`.
"""
function setup_env(proj::Project)
	empty!(Base.LOAD_PATH)
	push!(Base.LOAD_PATH, "@", "@stdlib")
	Pkg.activate(proj)
	return proj
end
setup_env(proj::ActiveProject) = setup_env(Project(proj))

#Last line

max_stars_repo_name: ma-laforge/JuliaApplications

max_issues_repo_name: ma-laforge/JuliaApplications

max_forks_repo_name: ma-laforge/JuliaApplications


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(cos(a+b*x), x) == :(b^-1*sin(a+b*x))
@test integrate(cos(a+b*x)^2, x) == :((1/2)*x+(1/2)*b^-1*cos(a+b*x)*sin(a+b*x))
@test integrate(cos(a+b*x)^3, x) == :(b^-1*sin(a+b*x)+-1//3*b^-1*sin(a+b*x)^3)
@test integrate(cos(a+b*x)^4, x) == :(3//8*x+1//4*b^-1*cos(a+b*x)^3*sin(a+b*x)+3//8*b^-1*cos(a+b*x)*sin(a+b*x))
@test integrate(cos(a+b*x)^5, x) == :(b^-1*sin(a+b*x)+-2//3*b^-1*sin(a+b*x)^3+1//5*b^-1*sin(a+b*x)^5)
@test integrate(cos(a+b*x)^6, x) == :(5//16*x+1//6*b^-1*cos(a+b*x)^5*sin(a+b*x)+5//16*b^-1*cos(a+b*x)*sin(a+b*x)+5//24*b^-1*cos(a+b*x)^3*sin(a+b*x))
@test integrate(cos(a+b*x)^7, x) == :(b^-1*sin(a+b*x)+-1*b^-1*sin(a+b*x)^3+-1//7*b^-1*sin(a+b*x)^7+3//5*b^-1*sin(a+b*x)^5)
@test integrate(cos(a+b*x)^8, x) == :(35//128*x+1//8*b^-1*cos(a+b*x)^7*sin(a+b*x)+7//48*b^-1*cos(a+b*x)^5*sin(a+b*x)+35//128*b^-1*cos(a+b*x)*sin(a+b*x)+35//192*b^-1*cos(a+b*x)^3*sin(a+b*x))
@test integrate(cos(a+b*x)^7//2, x) == :(10//21*b^-1*Elliptic.F((1/2)*a+(1/2)*b*x,2)+2//7*b^-1*cos(a+b*x)^5//2*sin(a+b*x)+10//21*b^-1*cos(a+b*x)^(1/2)*sin(a+b*x))
@test integrate(cos(a+b*x)^5//2, x) == :(6//5*b^-1*Elliptic.E((1/2)*a+(1/2)*b*x,2)+2//5*b^-1*cos(a+b*x)^3//2*sin(a+b*x))
@test integrate(cos(a+b*x)^3//2, x) == :(2//3*b^-1*Elliptic.F((1/2)*a+(1/2)*b*x,2)+2//3*b^-1*cos(a+b*x)^(1/2)*sin(a+b*x))
@test integrate(cos(a+b*x)^(1/2), x) == :(2*b^-1*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate(cos(a+b*x)^-1//2, x) == :(2*b^-1*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate(cos(a+b*x)^-3//2, x) == :(-2*b^-1*Elliptic.E((1/2)*a+(1/2)*b*x,2)+2*b^-1*cos(a+b*x)^-1//2*sin(a+b*x))
@test integrate(cos(a+b*x)^-5//2, x) == :(2//3*b^-1*Elliptic.F((1/2)*a+(1/2)*b*x,2)+2//3*b^-1*cos(a+b*x)^-3//2*sin(a+b*x))
@test integrate(cos(a+b*x)^-7//2, x) == :(-6//5*b^-1*Elliptic.E((1/2)*a+(1/2)*b*x,2)+2//5*b^-1*cos(a+b*x)^-5//2*sin(a+b*x)+6//5*b^-1*cos(a+b*x)^-1//2*sin(a+b*x))
@test integrate((c*cos(a+b*x))^7//2, x) == :(2//7*c*b^-1*(c*cos(a+b*x))^5//2*sin(a+b*x)+10//21*b^-1*c^3*(c*cos(a+b*x))^(1/2)*sin(a+b*x)+10//21*b^-1*c^4*(c*cos(a+b*x))^-1//2*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^5//2, x) == :(2//5*c*b^-1*(c*cos(a+b*x))^3//2*sin(a+b*x)+6//5*b^-1*c^2*(c*cos(a+b*x))^(1/2)*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^3//2, x) == :(2//3*c*b^-1*(c*cos(a+b*x))^(1/2)*sin(a+b*x)+2//3*b^-1*c^2*(c*cos(a+b*x))^-1//2*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^(1/2), x) == :(2*b^-1*(c*cos(a+b*x))^(1/2)*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^-1//2, x) == :(2*b^-1*(c*cos(a+b*x))^-1//2*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^-3//2, x) == :(2*b^-1*c^-1*(c*cos(a+b*x))^-1//2*sin(a+b*x)+-2*b^-1*c^-2*(c*cos(a+b*x))^(1/2)*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^-5//2, x) == :(2//3*b^-1*c^-1*(c*cos(a+b*x))^-3//2*sin(a+b*x)+2//3*b^-1*c^-2*(c*cos(a+b*x))^-1//2*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*cos(a+b*x))^-7//2, x) == :(2//5*b^-1*c^-1*(c*cos(a+b*x))^-5//2*sin(a+b*x)+6//5*b^-1*c^-3*(c*cos(a+b*x))^-1//2*sin(a+b*x)+-6//5*b^-1*c^-4*(c*cos(a+b*x))^(1/2)*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((a*cos(x)^2)^5//2, x) == :(1//5*(a*cos(x)^2)^5//2*tan(x)+4//15*a*(a*cos(x)^2)^3//2*tan(x)+8//15*a^2*(a*cos(x)^2)^(1/2)*tan(x))
@test integrate((a*cos(x)^2)^3//2, x) == :(1//3*(a*cos(x)^2)^3//2*tan(x)+2//3*a*(a*cos(x)^2)^(1/2)*tan(x))
@test integrate((a*cos(x)^2)^(1/2), x) == :((a*cos(x)^2)^(1/2)*tan(x))
@test integrate((a*cos(x)^2)^-1//2, x) == :((a*cos(x)^2)^-1//2*arctanh(sin(x))*cos(x))
@test integrate((a*cos(x)^2)^-3//2, x) == :((1/2)*a^-1*(a*cos(x)^2)^-1//2*tan(x)+(1/2)*a^-1*(a*cos(x)^2)^-1//2*arctanh(sin(x))*cos(x))
@test integrate((a*cos(x)^2)^-5//2, x) == :(1//4*a^-1*(a*cos(x)^2)^-3//2*tan(x)+3//8*a^-2*(a*cos(x)^2)^-1//2*tan(x)+3//8*a^-2*(a*cos(x)^2)^-1//2*arctanh(sin(x))*cos(x))
@test integrate((a*cos(x)^3)^5//2, x) == :(26//77*a^2*(a*cos(x)^3)^(1/2)*tan(x)+2//15*a^2*(a*cos(x)^3)^(1/2)*cos(x)^5*sin(x)+26//77*a^2*(a*cos(x)^3)^(1/2)*cos(x)^-3//2*Elliptic.F((1/2)*x,2)+26//165*a^2*(a*cos(x)^3)^(1/2)*cos(x)^3*sin(x)+78//385*a^2*(a*cos(x)^3)^(1/2)*cos(x)*sin(x))
@test integrate((a*cos(x)^3)^3//2, x) == :(14//45*a*(a*cos(x)^3)^(1/2)*sin(x)+2//9*a*(a*cos(x)^3)^(1/2)*cos(x)^2*sin(x)+14//15*a*(a*cos(x)^3)^(1/2)*cos(x)^-3//2*Elliptic.E((1/2)*x,2))
@test integrate((a*cos(x)^3)^(1/2), x) == :(2//3*(a*cos(x)^3)^(1/2)*tan(x)+2//3*(a*cos(x)^3)^(1/2)*cos(x)^-3//2*Elliptic.F((1/2)*x,2))
@test integrate((a*cos(x)^3)^-1//2, x) == :(-2*(a*cos(x)^3)^-1//2*cos(x)^3//2*Elliptic.E((1/2)*x,2)+2*(a*cos(x)^3)^-1//2*cos(x)*sin(x))
@test integrate((a*cos(x)^3)^-3//2, x) == :(10//21*a^-1*(a*cos(x)^3)^-1//2*sin(x)+2//7*a^-1*(a*cos(x)^3)^-1//2*sec(x)*tan(x)+10//21*a^-1*(a*cos(x)^3)^-1//2*cos(x)^3//2*Elliptic.F((1/2)*x,2))
@test integrate((a*cos(x)^3)^-5//2, x) == :(154//585*a^-2*(a*cos(x)^3)^-1//2*tan(x)+-154//195*a^-2*(a*cos(x)^3)^-1//2*cos(x)^3//2*Elliptic.E((1/2)*x,2)+2//13*a^-2*(a*cos(x)^3)^-1//2*sec(x)^4*tan(x)+22//117*a^-2*(a*cos(x)^3)^-1//2*sec(x)^2*tan(x)+154//195*a^-2*(a*cos(x)^3)^-1//2*cos(x)*sin(x))
@test integrate((a*cos(x)^4)^5//2, x) == :(63//256*a^2*(a*cos(x)^4)^(1/2)*tan(x)+1//10*a^2*(a*cos(x)^4)^(1/2)*cos(x)^7*sin(x)+9//80*a^2*(a*cos(x)^4)^(1/2)*cos(x)^5*sin(x)+21//128*a^2*(a*cos(x)^4)^(1/2)*cos(x)*sin(x)+21//160*a^2*(a*cos(x)^4)^(1/2)*cos(x)^3*sin(x)+63//256*x*a^2*(a*cos(x)^4)^(1/2)*sec(x)^2)
@test integrate((a*cos(x)^4)^3//2, x) == :(5//16*a*(a*cos(x)^4)^(1/2)*tan(x)+1//6*a*(a*cos(x)^4)^(1/2)*cos(x)^3*sin(x)+5//16*a*x*(a*cos(x)^4)^(1/2)*sec(x)^2+5//24*a*(a*cos(x)^4)^(1/2)*cos(x)*sin(x))
@test integrate((a*cos(x)^4)^(1/2), x) == :((1/2)*(a*cos(x)^4)^(1/2)*tan(x)+(1/2)*x*(a*cos(x)^4)^(1/2)*sec(x)^2)
@test integrate((a*cos(x)^4)^-1//2, x) == :((a*cos(x)^4)^-1//2*cos(x)*sin(x))
@test integrate((a*cos(x)^4)^-3//2, x) == :(a^-1*(a*cos(x)^4)^-1//2*cos(x)*sin(x)+1//5*a^-1*(a*cos(x)^4)^-1//2*sin(x)^2*tan(x)^3+2//3*a^-1*(a*cos(x)^4)^-1//2*sin(x)^2*tan(x))
@test integrate((a*cos(x)^4)^-5//2, x) == :(a^-2*(a*cos(x)^4)^-1//2*cos(x)*sin(x)+1//9*a^-2*(a*cos(x)^4)^-1//2*sin(x)^2*tan(x)^7+4//3*a^-2*(a*cos(x)^4)^-1//2*sin(x)^2*tan(x)+4//7*a^-2*(a*cos(x)^4)^-1//2*sin(x)^2*tan(x)^5+6//5*a^-2*(a*cos(x)^4)^-1//2*sin(x)^2*tan(x)^3)
@test integrate((c*cos(a+b*x)^m)^(m^-1), x) == :(b^-1*(c*cos(a+b*x)^m)^(m^-1)*tan(a+b*x))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^5, x) == :(30//77*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//11*b^-4*d^-1*(b*cos(c+d*x))^9//2*sin(c+d*x)+18//77*b^-2*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+30//77*b*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^4, x) == :(2//9*b^-3*d^-1*(b*cos(c+d*x))^7//2*sin(c+d*x)+14//15*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+14//45*b^-1*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^3, x) == :(10//21*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//7*b^-2*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^2, x) == :(2//5*b^-1*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+6//5*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x), x) == :(2//3*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//3*b*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2), x) == :(2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*sec(c+d*x), x) == :(2*b*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*sec(c+d*x)^2, x) == :(-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2*b*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*sec(c+d*x)^3, x) == :(2//3*b^2*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//3*b*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*sec(c+d*x)^4, x) == :(-6//5*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//5*b^3*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+6//5*b*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*sec(c+d*x)^5, x) == :(2//7*b^4*d^-1*(b*cos(c+d*x))^-7//2*sin(c+d*x)+10//21*b^2*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+10//21*b*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*sec(c+d*x)^6, x) == :(-14//15*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//9*b^5*d^-1*(b*cos(c+d*x))^-9//2*sin(c+d*x)+14//15*b*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+14//45*b^3*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^4, x) == :(2//11*b^-3*d^-1*(b*cos(c+d*x))^9//2*sin(c+d*x)+18//77*b^-1*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+30//77*b*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+30//77*b^2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^3, x) == :(14//45*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+2//9*b^-2*d^-1*(b*cos(c+d*x))^7//2*sin(c+d*x)+14//15*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^2, x) == :(2//7*b^-1*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+10//21*b^2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x), x) == :(2//5*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+6//5*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2, x) == :(2//3*b*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//3*b^2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x), x) == :(2*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x)^2, x) == :(2*b^2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x)^3, x) == :(2*b^2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-2*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x)^4, x) == :(2//3*b^3*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//3*b^2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x)^5, x) == :(2//5*b^4*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+6//5*b^2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-6//5*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x)^6, x) == :(2//7*b^5*d^-1*(b*cos(c+d*x))^-7//2*sin(c+d*x)+10//21*b^3*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+10//21*b^2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^3//2*sec(c+d*x)^7, x) == :(2//9*b^6*d^-1*(b*cos(c+d*x))^-9//2*sin(c+d*x)+14//15*b^2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+14//45*b^4*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+-14//15*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^3, x) == :(18//77*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+2//11*b^-2*d^-1*(b*cos(c+d*x))^9//2*sin(c+d*x)+30//77*b^2*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+30//77*b^3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^2, x) == :(2//9*b^-1*d^-1*(b*cos(c+d*x))^7//2*sin(c+d*x)+14//45*b*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+14//15*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x), x) == :(2//7*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b^2*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+10//21*b^3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2, x) == :(2//5*b*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+6//5*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x), x) == :(2//3*b^2*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//3*b^3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^2, x) == :(2*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^3, x) == :(2*b^3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^4, x) == :(2*b^3*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-2*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^5, x) == :(2//3*b^4*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//3*b^3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^6, x) == :(2//5*b^5*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+6//5*b^3*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-6//5*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^7, x) == :(2//7*b^6*d^-1*(b*cos(c+d*x))^-7//2*sin(c+d*x)+10//21*b^4*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+10//21*b^3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^5//2*sec(c+d*x)^8, x) == :(2//9*b^7*d^-1*(b*cos(c+d*x))^-9//2*sin(c+d*x)+14//15*b^3*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+14//45*b^5*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+-14//15*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^7//2, x) == :(2//7*b*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b^3*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+10//21*b^4*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^6, x) == :(2//11*b^-5*d^-1*(b*cos(c+d*x))^9//2*sin(c+d*x)+18//77*b^-3*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+30//77*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+30//77*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^5, x) == :(2//9*b^-4*d^-1*(b*cos(c+d*x))^7//2*sin(c+d*x)+14//45*b^-2*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+14//15*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^4, x) == :(2//7*b^-3*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+10//21*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^3, x) == :(2//5*b^-2*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+6//5*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^2, x) == :(2//3*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x), x) == :(2*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2, x) == :(2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*sec(c+d*x), x) == :(2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-2*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*sec(c+d*x)^2, x) == :(2//3*b*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*sec(c+d*x)^3, x) == :(6//5*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+2//5*b^2*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+-6//5*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*sec(c+d*x)^4, x) == :(2//7*b^3*d^-1*(b*cos(c+d*x))^-7//2*sin(c+d*x)+10//21*b*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+10//21*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-1//2*sec(c+d*x)^5, x) == :(14//15*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+2//9*b^4*d^-1*(b*cos(c+d*x))^-9//2*sin(c+d*x)+14//45*b^2*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+-14//15*b^-1*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^7, x) == :(2//11*b^-6*d^-1*(b*cos(c+d*x))^9//2*sin(c+d*x)+18//77*b^-4*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+30//77*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+30//77*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^6, x) == :(2//9*b^-5*d^-1*(b*cos(c+d*x))^7//2*sin(c+d*x)+14//45*b^-3*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+14//15*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^5, x) == :(2//7*b^-4*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+10//21*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^4, x) == :(2//5*b^-3*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+6//5*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^3, x) == :(2//3*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//3*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^2, x) == :(2*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x), x) == :(2*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2, x) == :(2*b^-1*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-2*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*sec(c+d*x), x) == :(2//3*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//3*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*sec(c+d*x)^2, x) == :(2//5*b*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+6//5*b^-1*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-6//5*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*sec(c+d*x)^3, x) == :(10//21*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//7*b^2*d^-1*(b*cos(c+d*x))^-7//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-3//2*sec(c+d*x)^4, x) == :(2//9*b^3*d^-1*(b*cos(c+d*x))^-9//2*sin(c+d*x)+14//15*b^-1*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+14//45*b*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+-14//15*b^-2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^8, x) == :(2//11*b^-7*d^-1*(b*cos(c+d*x))^9//2*sin(c+d*x)+18//77*b^-5*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+30//77*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+30//77*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^7, x) == :(2//9*b^-6*d^-1*(b*cos(c+d*x))^7//2*sin(c+d*x)+14//45*b^-4*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+14//15*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^6, x) == :(2//7*b^-5*d^-1*(b*cos(c+d*x))^5//2*sin(c+d*x)+10//21*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+10//21*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^5, x) == :(2//5*b^-4*d^-1*(b*cos(c+d*x))^3//2*sin(c+d*x)+6//5*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^4, x) == :(2//3*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*sin(c+d*x)+2//3*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^3, x) == :(2*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^2, x) == :(2*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x), x) == :(2*b^-2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-2*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2, x) == :(2//3*b^-1*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+2//3*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*sec(c+d*x), x) == :(2//5*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+6//5*b^-2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-6//5*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*sec(c+d*x)^2, x) == :(2//7*b*d^-1*(b*cos(c+d*x))^-7//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*cos(c+d*x))^-3//2*sin(c+d*x)+10//21*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-5//2*sec(c+d*x)^3, x) == :(14//45*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+2//9*b^2*d^-1*(b*cos(c+d*x))^-9//2*sin(c+d*x)+14//15*b^-2*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-14//15*b^-3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^-7//2, x) == :(2//5*b^-1*d^-1*(b*cos(c+d*x))^-5//2*sin(c+d*x)+6//5*b^-3*d^-1*(b*cos(c+d*x))^-1//2*sin(c+d*x)+-6//5*b^-4*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^7//2, x) == :(3//8*x*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2+1//4*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^5//2*sin(c+d*x)+3//8*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^5//2, x) == :(d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)+-1//3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^3//2, x) == :((1/2)*x*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2+(1/2)*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2), x) == :(d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2, x) == :(x*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2)
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2, x) == :(d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2, x) == :(d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^-7//2, x) == :((1/2)*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2*sin(c+d*x)+(1/2)*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^-9//2, x) == :(d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2*sin(c+d*x)+1//3*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-7//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^(1/2)*cos(c+d*x)^-11//2, x) == :(1//4*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-9//2*sin(c+d*x)+3//8*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2*sin(c+d*x)+3//8*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^5//2, x) == :(3//8*b*x*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2+1//4*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^5//2*sin(c+d*x)+3//8*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^3//2, x) == :(b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)+-1//3*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^(1/2), x) == :((1/2)*b*x*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2+(1/2)*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-1//2, x) == :(b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-3//2, x) == :(b*x*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2)
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-5//2, x) == :(b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-7//2, x) == :(b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-9//2, x) == :((1/2)*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2*sin(c+d*x)+(1/2)*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-11//2, x) == :(b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2*sin(c+d*x)+1//3*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-7//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^3//2*cos(c+d*x)^-13//2, x) == :(1//4*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-9//2*sin(c+d*x)+3//8*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2*sin(c+d*x)+3//8*b*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^5//2, x) == :(b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)+-2//3*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)^3+1//5*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)^5)
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^3//2, x) == :(3//8*x*b^2*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2+1//4*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^5//2*sin(c+d*x)+3//8*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^(1/2), x) == :(b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)+-1//3*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-1//2, x) == :((1/2)*x*b^2*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2+(1/2)*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-3//2, x) == :(b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-5//2, x) == :(x*b^2*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2)
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-7//2, x) == :(b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-9//2, x) == :(b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-11//2, x) == :((1/2)*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2*sin(c+d*x)+(1/2)*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-13//2, x) == :(b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-3//2*sin(c+d*x)+1//3*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-7//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^5//2*cos(c+d*x)^-15//2, x) == :(1//4*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-9//2*sin(c+d*x)+3//8*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-5//2*sin(c+d*x)+3//8*b^2*d^-1*(b*cos(c+d*x))^(1/2)*cos(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^11//2, x) == :(d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)+-2//3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)^3+1//5*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)^5)
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^9//2, x) == :(3//8*x*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)+1//4*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^7//2*sin(c+d*x)+3//8*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^7//2, x) == :(d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)+-1//3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^5//2, x) == :((1/2)*x*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)+(1/2)*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2, x) == :(d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2), x) == :(x*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2, x) == :(d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2, x) == :(d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^-5//2, x) == :((1/2)*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x))+(1/2)*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^-7//2, x) == :(d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2*sin(c+d*x)+1//3*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-5//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^-1//2*cos(c+d*x)^-9//2, x) == :(1//4*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-7//2*sin(c+d*x)+3//8*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x))+3//8*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^11//2, x) == :(3//8*x*b^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)+1//4*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^7//2*sin(c+d*x)+3//8*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^9//2, x) == :(b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)+-1//3*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^7//2, x) == :((1/2)*x*b^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)+(1/2)*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^5//2, x) == :(b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^3//2, x) == :(x*b^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^(1/2), x) == :(b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^-1//2, x) == :(b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^-3//2, x) == :((1/2)*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x))+(1/2)*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^-5//2, x) == :(b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2*sin(c+d*x)+1//3*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-5//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^-3//2*cos(c+d*x)^-7//2, x) == :(1//4*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-7//2*sin(c+d*x)+3//8*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x))+3//8*b^-1*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^13//2, x) == :(3//8*x*b^-2*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)+1//4*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^7//2*sin(c+d*x)+3//8*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^11//2, x) == :(b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)+-1//3*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^9//2, x) == :((1/2)*x*b^-2*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)+(1/2)*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^7//2, x) == :(b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^5//2, x) == :(x*b^-2*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^3//2, x) == :(b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x)))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^(1/2), x) == :(b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^-1//2, x) == :((1/2)*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x))+(1/2)*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^-3//2, x) == :(b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-1//2*sin(c+d*x)+1//3*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-5//2*sin(c+d*x)^3)
@test integrate((b*cos(c+d*x))^-5//2*cos(c+d*x)^-5//2, x) == :(1//4*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-7//2*sin(c+d*x)+3//8*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^(1/2)*arctanh(sin(c+d*x))+3//8*b^-2*d^-1*(b*cos(c+d*x))^-1//2*cos(c+d*x)^-3//2*sin(c+d*x))
@test integrate(csc(a+b*x)^(1/2)*cos(a+b*x), x) == :(2*b^-1*csc(a+b*x)^-1//2)
@test integrate(csc(a+b*x)^-1//2*cos(a+b*x), x) == :(2//3*b^-1*csc(a+b*x)^-3//2)
@test integrate(cos(a+b*x)^2*csc(a+b*x)^(1/2), x) == :(2//3*b^-1*csc(a+b*x)^-1//2*cos(a+b*x)+4//3*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(cos(a+b*x)^2*csc(a+b*x)^-1//2, x) == :(2//5*b^-1*csc(a+b*x)^-3//2*cos(a+b*x)+4//5*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(cos(x)^3*csc(x)^9//2, x) == :(-2//7*csc(x)^7//2+2//3*csc(x)^3//2)
@test integrate(cos(a+b*x)^3*csc(a+b*x)^(1/2), x) == :(2*b^-1*csc(a+b*x)^-1//2+-2//5*b^-1*csc(a+b*x)^-5//2)
@test integrate(cos(a+b*x)^3*csc(a+b*x)^-1//2, x) == :(-2//7*b^-1*csc(a+b*x)^-7//2+2//3*b^-1*csc(a+b*x)^-3//2)
@test integrate(cos(a+b*x)^4*csc(a+b*x)^(1/2), x) == :(2//7*b^-1*cos(a+b*x)^3*csc(a+b*x)^-1//2+4//7*b^-1*csc(a+b*x)^-1//2*cos(a+b*x)+8//7*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(cos(a+b*x)^4*csc(a+b*x)^-1//2, x) == :(2//9*b^-1*cos(a+b*x)^3*csc(a+b*x)^-3//2+4//15*b^-1*csc(a+b*x)^-3//2*cos(a+b*x)+8//15*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(x)^7//3*cos(x), x) == :(-3//4*csc(x)^4//3)
@test integrate(csc(a+b*x)^(1/2)*sec(a+b*x), x) == :(b^-1*arctanh(csc(a+b*x)^(1/2))+-1*b^-1*arctan(csc(a+b*x)^(1/2)))
@test integrate(csc(a+b*x)^-1//2*sec(a+b*x), x) == :(b^-1*arctan(csc(a+b*x)^(1/2))+b^-1*arctanh(csc(a+b*x)^(1/2)))
@test integrate(csc(a+b*x)^(1/2)*sec(a+b*x)^2, x) == :(b^-1*csc(a+b*x)^-1//2*sec(a+b*x)+b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^-1//2*sec(a+b*x)^2, x) == :(b^-1*csc(a+b*x)^-3//2*sec(a+b*x)+-1*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^(1/2)*sec(a+b*x)^3, x) == :(-3//4*b^-1*arctan(csc(a+b*x)^(1/2))+3//4*b^-1*arctanh(csc(a+b*x)^(1/2))+(1/2)*b^-1*csc(a+b*x)^-1//2*sec(a+b*x)^2)
@test integrate(csc(a+b*x)^-1//2*sec(a+b*x)^3, x) == :(1//4*b^-1*arctan(csc(a+b*x)^(1/2))+1//4*b^-1*arctanh(csc(a+b*x)^(1/2))+(1/2)*b^-1*csc(a+b*x)^-3//2*sec(a+b*x)^2)
@test integrate(csc(a+b*x)^(1/2)*sec(a+b*x)^4, x) == :(1//3*b^-1*csc(a+b*x)^-1//2*sec(a+b*x)^3+5//6*b^-1*csc(a+b*x)^-1//2*sec(a+b*x)+5//6*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*x,2))
@test integrate(csc(a+b*x)^-1//2*sec(a+b*x)^4, x) == :((1/2)*b^-1*csc(a+b*x)^-3//2*sec(a+b*x)+1//3*b^-1*csc(a+b*x)^-3//2*sec(a+b*x)^3+-1//2*b^-1*csc(a+b*x)^(1/2)*sin(a+b*x)^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*x,2))

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.0 (a cos)^m (b trg)^n.jl


text:
# Vendoring `Base.afoldl`
afoldl(op,a) = a
afoldl(op,a,b) = op(a,b)
afoldl(op,a,b,c...) = afoldl(op, op(a,b), c...)
function afoldl(op,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,qs...)
    y = op(op(op(op(op(op(op(op(op(op(op(op(op(op(op(a,b),c),d),e),f),g),h),i),j),k),l),m),n),o),p)
    for x in qs; y = op(y,x); end
    y
end

@inline getindex(x, inds...) = Base.getindex(x, inds...)
@def getindex(xs::Tuple, ind::UnitRange) = ntuple(i -> xs[ind[i]], length(ind))

@inline filter(f, xs) = Base.filter(f, xs)
@def filter(f::F, xs::Tuple) where {F} =
    afoldl((ys, x) -> f(x) ? (ys..., x) : ys, (), xs...)

zip(x, xs...) = Base.zip(x, xs...)

# Not using `zip(tuples::NTuple{N,Any}...)` to avoid having unbound
# type parameter.
@def function zip(x::NTuple{N,Any}, xs::NTuple{N,Any}...) where {N}
    tuples = (x, xs...)
    return ntuple(i -> map(t -> t[i], tuples), N)
end

@inline flatten(itr) = Iterators.flatten(itr)

@def flatten(xs::Tuple) = _flatten(xs...)
@inline _flatten() = ()
@inline _flatten(x::Tuple, xs::Tuple...) = (x..., _flatten(xs...)...)

@inline enumerate(itr) = Base.enumerate(itr)

@def enumerate(xs::Tuple) = Specialized.zip(ntuple(identity, length(xs)), xs)

struct _InitialValue end

struct BottomRF{T}
    rf::T
end

@inline (op::BottomRF)(::_InitialValue, x) = Base.reduce_first(op.rf, x)
@inline (op::BottomRF)(acc, x) = op.rf(acc, x)

@inline accumulate(args...; kw...) = Base.accumulate(args...; kw...)

@def function accumulate(op::F, xs::Tuple; init = _InitialValue()) where {F}
    rf = BottomRF(op)
    ys, = afoldl(((), init), xs...) do (ys, acc), x
        acc = rf(acc, x)
        (ys..., acc), acc
    end
    return ys
end

@inline cumsum(args...; kw...) = Base.cumsum(args...; kw...)
@def cumsum(xs::Tuple) = Specialized.accumulate(Base.add_sum, xs)

@inline cumprod(args...; kw...) = Base.cumprod(args...; kw...)
@def cumprod(xs::Tuple) = Specialized.accumulate(Base.mul_prod, xs)

@inline foreach(f, xs, rest...) = Base.foreach(f, xs, rest...)

# Not using `@def` so that `Specialized.foreach` is always defined via `foldl`.
@inline foreach(f, xs::Tuple, rest::Tuple...) = Specialized.foreach(f, xs, rest...)
@inline Specialized.foreach(f, xs) = afoldl((_, x) -> (f(x); nothing), nothing, xs...)
@inline Specialized.foreach(f, xs, rest...) =
    afoldl((_, args) -> (f(args...); nothing), nothing, zip(xs, rest...)...)
# https://github.com/JuliaLang/julia/pull/31901

@inline any(f, itr) = Base.any(f, itr)
@inline any(itr) = Base.any(itr)

@def any(xs::Tuple) = _any(identity, xs...)
@def any(f::F, xs::Tuple) where {F} = _any(f, xs...)
@inline _any(f) = false
@inline function _any(f::F, x, xs...) where {F}
    y = f(x)
    y === true && return true
    return y | _any(f, xs...)
end

@inline all(f, itr) = Base.all(f, itr)
@inline all(itr) = Base.all(itr)

@def all(xs::Tuple) = _all(identity, xs...)
@def all(f::F, xs::Tuple) where {F} = _all(f, xs...)
@inline _all(f) = true
@inline function _all(f::F, x, xs...) where {F}
    y = f(x)
    y === false && return false
    return y & _all(f, xs...)
end

@inline findfirst(xs) = Base.findfirst(xs)
@inline findfirst(f, xs) = Base.findfirst(f, xs)
@def findfirst(xs::Tuple{Vararg{Bool}}) = _findfirst(identity, 1, xs...)
@def findfirst(f::F, xs::Tuple) where {F} = _findfirst(f, 1, xs...)
@inline _findfirst(f, n) = nothing
@inline _findfirst(f::F, n, x, xs...) where {F} = f(x) ? n : _findfirst(f, n + 1, xs...)

@inline findlast(xs) = Base.findlast(xs)
@inline findlast(f, xs) = Base.findlast(f, xs)
@def findlast(xs::Tuple{Vararg{Bool}}) = _findlast(identity, 1, xs...)
@def findlast(f::F, xs::Tuple) where {F} = _findlast(f, 1, xs...)
@inline _findlast(f, n) = nothing
@inline function _findlast(f::F, n, x, xs...) where {F}
    i = _findlast(f, n + 1, xs...)
    i === nothing || return i
    return f(x) ? n : nothing
end

@inline findall(xs) = Base.findall(xs)
@inline findall(f, xs) = Base.findall(f, xs)
@def findall(xs::Tuple{Vararg{Bool}}) = _findall(identity, 1, xs...)
@def findall(f::F, xs::Tuple) where {F} = _findall(f, 1, xs...)
@inline _findall(f, n) = ()
@inline _findall(f, n, x, xs...) = ((f(x) ? (n,) : ())..., _findall(f, n + 1, xs...)...)

@inline findmax(args...; kw...) = Base.findmax(args...; kw...)
@def findmax(xs::Tuple{}) = throw(ArgumentError("collection must be non-empty"))
@def findmax(xs::Tuple) = _findmax(zip(xs, ntuple(identity, length(xs)))...)
@inline _findmax(ans) = ans
@inline function _findmax((x, i), (y, j), rest...)
    x != x && return (x, i)
    return _findmax(y != y || isless(x, y) ? (y, j) : (x, i), rest...)
end

@inline findmin(args...; kw...) = Base.findmin(args...; kw...)
@def findmin(xs::Tuple{}) = throw(ArgumentError("collection must be non-empty"))
@def findmin(xs::Tuple) = _findmin(zip(xs, ntuple(identity, length(xs)))...)
@inline _findmin(ans) = ans
@inline function _findmin((x, i), (y, j), rest...)
    x != x && return (x, i)
    return _findmin(y != y || isless(y, x) ? (y, j) : (x, i), rest...)
end

@inline argmax(args...; kw...) = Base.argmax(args...; kw...)
@def argmax(xs::Tuple) = Specialized.findmax(xs)[2]

@inline argmin(args...; kw...) = Base.argmin(args...; kw...)
@def argmin(xs::Tuple) = Specialized.findmin(xs)[2]

@inline maximum(args...; kw...) = Base.maximum(args...; kw...)
@def maximum(xs::Tuple) = Specialized.findmax(xs)[1]

@inline minimum(args...; kw...) = Base.minimum(args...; kw...)
@def minimum(xs::Tuple) = Specialized.findmin(xs)[1]

@inline extrema(args...; kw...) = Base.extrema(args...; kw...)
@def extrema(xs::Tuple{}) = throw(ArgumentError("collection must be non-empty"))
@def extrema(xs::Tuple) = Specialized.extrema(identity, xs)
@def extrema(f, xs::Tuple{}) = throw(ArgumentError("collection must be non-empty"))
@def function extrema(f::F, xs::Tuple) where {F}
    x = xs[1]
    x != x && return (x, x)
    return _extrema(f, (x, x), Base.tail(xs)...)
end
@inline _extrema(f, ans) = ans
@inline function _extrema(f::F, (l, g), x, xs...) where {F}
    x != x && return (x, x)
    return _extrema(f, (isless(x, l) ? x : l, isless(g, x) ? x : g), xs...)
end

@inline in(x, xs) = Base.in(x, xs)
@inline in(x) = Base.Fix2(in, x)
@inline Specialized.in(x) = Base.Fix2(Specialized.in, x)
@def in(x, xs::Tuple) = _any(==(x), xs...)

@inline unique(f, itr) = Base.unique(f, itr)
@inline unique(itr) = Base.unique(itr)
@def unique(xs::Tuple{}) = ()
@def unique(xs::Tuple) =
    afoldl((), xs...) do xs′, x
        x in xs′ ? xs′ : (xs′..., x)
    end
@def unique(f, xs::Tuple{}) = ()
@def function unique(f::F, xs::Tuple) where {F}
    ys = map(f, xs)
    xs′, = afoldl(((), ()), zip(xs, ys)...) do (xs, ys), (x, y)
        y in ys ? (xs, ys) : ((ys..., y), (xs..., x))
    end
    return xs′
end

@inline union(args...) = Base.union(args...)

@def union(a::Tuple) = Specialized.unique(a)
@def union(a::Tuple, b::Tuple) = unique((a..., b...))
@def union(a::Tuple, bs::Tuple...) = afoldl(union, a, bs...)

@inline intersect(s, itrs...) = Base.intersect(s, itrs...)
@def intersect(a::Tuple) = a
@def intersect(a::Tuple{}, b::Tuple{}) = ()
@def intersect(a::Tuple, b::Tuple{}) = ()
@def intersect(a::Tuple{}, b::Tuple) = ()
@def intersect(a::Tuple, b::Tuple) = unique(filter(x -> x in b, a))

@def intersect(a::Tuple, b::Tuple, cs::Tuple...) =
    foldl(Specialized.intersect, cs, init = Specialized.intersect(a, b))

@inline setdiff(args...) = Base.setdiff(args...)
@def setdiff(a::Tuple, b::Tuple{}) = a
@def setdiff(a::Tuple, b::Tuple) = setdiff(_exclude(a, b[1]), Base.tail(b))
@inline _exclude(a, b) = Base.foldl((ys, x) -> x == b ? ys : (ys..., x), a; init = ())

@def setdiff(a::Tuple, b::Tuple, cs::Tuple...) =
    Base.foldl(Specialized.setdiff, cs, init = Specialized.setdiff(a, b))

@def symdiff(a::Tuple, b::Tuple) =
    foldl(b; init = a) do a, x
        i = findfirst(==(x), a)
        i == nothing ? (a..., x) : _deleteat(a, i)
    end
@def symdiff(a::Tuple, b::Tuple, cs::Tuple...) =
    foldl(Specialized.symdiff, cs, init = Specialized.symdiff(a, b))

_deleteat(a, i::Int) =
    foldl(a; init = ((), 1)) do (ys, j), x
        ((i == j ? ys : (ys..., x)), j + 1)
    end

@inline issubset(a, b) = Base.issubset(a, b)
@def issubset(a::Tuple, b::Tuple) = Specialized.all(x -> x in b, a)

if VERSION >= v"1.5"
    @inline isdisjoint(a, b) = Base.isdisjoint(a, b)
else
    @inline isdisjoint(a, b) = isempty(intersect(a, b))
end
@def isdisjoint(a::Tuple, b::Tuple) = !Specialized.any(x -> x in b, a)

sort(v; kw...) = Base.sort(v; kw...)

@def sort(
    v::Tuple;
    lt = isless,
    by = identity,
    rev::Union{Bool,Nothing} = nothing,
    order = Base.Forward,
) = _sort(Base.ord(lt, by, rev, order), v)

@inline _sort(order, ::Tuple{}) = ()
@inline _sort(order, x::Tuple{Any}) = x
@inline _sort(order, (x, y)::Tuple{Any,Any}) = Base.lt(order, y, x) ? (y, x) : (x, y)
@inline function _sort(order, v)
    left, right = _halve(v)
    return _mergesorted(order, _sort(order, left), _sort(order, right))
end

@inline _mergesorted(order, ::Tuple{}, ::Tuple{}) = ()
@inline _mergesorted(order, ::Tuple{}, right) = right
@inline _mergesorted(order, left, ::Tuple{}) = left
@inline function _mergesorted(order, left, right)
    a = left[1]
    b = right[1]
    if Base.lt(order, b, a)
        return (b, _mergesorted(order, left, Base.tail(right))...)
    else
        return (a, _mergesorted(order, Base.tail(left), right)...)
    end
end

@inline function _halve(v::NTuple{N,Any}) where {N}
    m = N ÷ 2
    return (getindex(v, 1:m), getindex(v, m+1:N))
end

# Compilation takes too long (> 10 sec) for `length(v::Tuple) > 13`.
const Any14{N} =
    Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,Any,Any,Any,Any,Any,Vararg{Any,N}}
@nospecialize
sort(v::Any14; kw...) = Tuple(Base.sort!(collect(v); kw...))
@specialize


text:
@json struct BranchChecks
    enforcement_level::String
    contexts::Vector{String}
end

@json struct BranchProtection
    enabled::Bool
    required_status_checks::BranchChecks
end

@json struct BranchLinks
    html::String
    self::String
end

@json struct Branch
    name::String
    commit::Commit
    _links => links::BranchLinks
    protected::Bool
    protection::BranchProtection
    protection_url::String
end

endpoint(::GitHubAPI, ::typeof(get_branch), owner::AStr, repo::AStr, branch::AStr) =
    Endpoint(:GET, "/repos/$owner/$repo/branches/$branch")
into(::GitHubAPI, ::typeof(get_branch)) = Branch

# Get Repository Branches
# https://docs.github.com/en/rest/reference/repos#list-branches
function endpoint(
    ::GitHubAPI,
    ::typeof(get_branches),
    owner::AStr,
    repo::AStr,
)
    return Endpoint(:GET, "repos/$owner/$repo/branches")
end
into(::GitHubAPI, ::typeof(get_branches)) = Vector{Branch}

# Delete Branch (Delete Reference)
# https://docs.github.com/en/rest/reference/git#delete-a-reference
# The :ref in the URL must be formatted as heads/<branch name> for branches
function endpoint(
    ::GitHubAPI,
    ::typeof(delete_branch),
    owner::AStr,
    repo::AStr,
    branch::AStr,
)
    return Endpoint(:DELETE, "/repos/$owner/$repo/git/refs/heads/$branch")
end
postprocessor(::GitHubAPI, ::typeof(delete_branch)) = DoNothing()


text:
export AsyncPipelineStage

import Rocket: async

"""
    AsyncPipelineStage <: AbstractPipelineStage

Applies the `async()` operator from `Rocket.jl` library to the given pipeline

See also: [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref), [`EmptyPipelineStage`](@ref), [`CompositePipelineStage`](@ref)
"""
struct AsyncPipelineStage <: AbstractPipelineStage end

apply_pipeline_stage(::AsyncPipelineStage, factornode, tag, stream) = stream |> async()


text:
#The scope of this test set is to check if all formulations provide more or less
#the same result #if errors are added on the measurements.

@testset "Measurements with errors - comparison between formulations" begin

    # set measurement path
    msr_path = joinpath(mktempdir(),"temp.csv")
    data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
    if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
    if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

    # insert the load profiles
    _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

    # transform data model
    data = _PMD.transform_data_model(data);

    @testset "ACP input-rwlav" begin
        # set model
        crit     = "rwlav"
        model    = _PMD.ACPUPowerModel #this is going to be used for the SE

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vi","vr"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
        _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        original_se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)
        delta_ref, max_ref, avg_ref = _PMDSE.calculate_voltage_magnitude_error(original_se_result, pf_result)
        for data_model in [_PMD.ACRUPowerModel,_PMD.IVRUPowerModel]

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, data_model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max-max_ref, 0.0; atol = 1e-3)
            @test isapprox(avg-avg_ref, 0.0; atol = 1e-3)
        end
    end
    @testset "ACP input-WLS" begin
        # set model
        crit     = "wls"
        model    = _PMD.ACPUPowerModel #this is going to be used for the SE

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vi","vr"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
        _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        original_se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)
        delta_ref, max_ref, avg_ref = _PMDSE.calculate_voltage_magnitude_error(original_se_result, pf_result)
        for data_model in [_PMD.ACRUPowerModel,_PMD.IVRUPowerModel]

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, data_model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max-max_ref, 0.0; atol = 2e-3)
            @test isapprox(avg-avg_ref, 0.0; atol = 1e-3)
        end
    end
    @testset "ACR input-rwlav" begin
        # set model
        crit     = "rwlav"
        model    = _PMD.ACRUPowerModel #this is going to be used for the SE

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vi","vr"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
        _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        original_se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)
        delta_ref, max_ref, avg_ref = _PMDSE.calculate_voltage_magnitude_error(original_se_result, pf_result)
        for data_model in [_PMD.ACPUPowerModel,_PMD.IVRUPowerModel]

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, data_model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max-max_ref, 0.0; atol = 1e-3)
            @test isapprox(avg-avg_ref, 0.0; atol = 1e-3)
        end
    end
    @testset "ACR input-WLS" begin
        # set model
        crit     = "wls"
        model    = _PMD.ACRUPowerModel #this is going to be used for the SE

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vi","vr"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
        _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 10)

        # solve the state estimation
        original_se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)
        delta_ref, max_ref, avg_ref = _PMDSE.calculate_voltage_magnitude_error(original_se_result, pf_result)
        for data_model in [_PMD.ACPUPowerModel,_PMD.IVRUPowerModel]

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, data_model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max-max_ref, 0.0; atol = 1e-3)
            @test isapprox(avg-avg_ref, 0.0; atol = 1e-3)
        end
    end
    @testset "IVR input-rwlav" begin
        # set model
        crit     = "rwlav"
        model    = _PMD.IVRUPowerModel #this is going to be used for the SE

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vi","vr"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
        _PMDSE.add_voltage_measurement!(data, pf_result, 0.0005)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        original_se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)
        delta_ref, max_ref, avg_ref = _PMDSE.calculate_voltage_magnitude_error(original_se_result, pf_result)
        for data_model in [_PMD.ACPUPowerModel,_PMD.ACRUPowerModel]

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, data_model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max-max_ref, 0.0; atol = 7e-3)
            @test isapprox(avg-avg_ref, 0.0; atol = 2.6e-3) #make 1.1!
        end
    end
    @testset "IVR input-WLS" begin
        # set model
        crit     = "wls"
        model    = _PMD.IVRUPowerModel #this is going to be used for the SE

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vi","vr"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
        _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 10)

        # solve the state estimation
        original_se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)
        delta_ref, max_ref, avg_ref = _PMDSE.calculate_voltage_magnitude_error(original_se_result, pf_result)
        for data_model in [_PMD.ACPUPowerModel,_PMD.ACRUPowerModel]

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, data_model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max-max_ref, 0.0; atol = 1e-3)
            @test isapprox(avg-avg_ref, 0.0; atol = 1e-3)
        end
    end
end


text:
module OpCodes

using OffsetArrays
import Base: iterate, setindex!, getindex

# process an instruction to an opcode
function process_instruction!(digit, instruction)
    digits!(digit, instruction)
    Val(instruction % 100)
end

# The IntComputer holding all the info
struct IntComputer
    data::OffsetArray{Int, 1}
    nargs::Dict{Type, Int}
    input::Channel{Int}
    output::Channel{Int}
    relativeindex::Ref{Int}
    extradata::Dict{Int, Int}
    finished::Ref{Bool}
    moreargs::Dict{Any, Any}
end
function IntComputer(data, input = Channel{Int}(Inf), output = Channel{Int}(Inf))
    nargs = Dict(x.sig.parameters[3] => x.nargs-3 for x in methods(run_opcode!).ms)
    IntComputer(OffsetArray(copy(data), 0:length(data)-1), nargs, input, output, Ref(0), Dict{Int, Int}(), Ref(false), Dict{Any, Any}())
end
IntComputer(str::AbstractString, input = Channel{Int}(Inf), output = Channel{Int}(Inf)) =
    IntComputer(parse.(Int, split(str, ',')), input, output)

# Handle getindex and setindex for different modes
struct ModeIndex{T}
    arg::Int
end
ModeIndex(T, arg) = ModeIndex{T}(arg)

setindex!(ic::IntComputer, val::Int, mi::ModeIndex) = setindex!(ic, val, mi.arg)
setindex!(ic::IntComputer, val::Int, mi::ModeIndex{2}) = setindex!(ic, val, ic.relativeindex[] + mi.arg)
setindex!(ic::IntComputer, val::Int, arg) = arg > length(ic.data)-1 ? ic.extradata[arg] = val : ic.data[arg] = val
getindex(ic::IntComputer, mi::ModeIndex{0}) = getindex(ic, mi.arg)
getindex(ic::IntComputer, mi::ModeIndex{1}) = mi.arg
getindex(ic::IntComputer, mi::ModeIndex{2}) = getindex(ic, ic.relativeindex[] + mi.arg)
getindex(ic::IntComputer, arg) = arg > length(ic.data)-1 ? get!(ic.extradata, arg, 0) : ic.data[arg]

addinput!(ic, inp) = put!(ic.input, inp)
takeinput!(ic) = take!(ic.input)
getoutput(ic) = (close(ic.output); collect(ic.output))

compute!(ic::IntComputer, input) = (addinput!(ic, input); compute!(ic))
function compute!(ic::IntComputer)
    state = 0
    ic.relativeindex[] = 0
    modeargs = Vector{ModeIndex}(undef, 3)
    modes = Vector{Int}(undef, 5)
    while state > -2
        opcode = process_instruction!(modes, ic.data[state])
        nargs = ic.nargs[typeof(opcode)]
        resize!(modeargs, nargs)
        modeargs .= ModeIndex.(modes[2 .+ (1:nargs)], ic.data[state .+ (1:nargs)])
        f = run_opcode!(ic, opcode, modeargs...)
        state = f == -1 ? state + nargs + 1 : f
    end
    ic.finished[] = true
    ic
end

run_opcode!(::Any, ::Val) = error("unknown opcode")
run_opcode!(ic, ::Val{99}) = -2
run_opcode!(ic, ::Val{1}, src1, src2, dest) = (ic[dest] = ic[src1] + ic[src2]; -1)
run_opcode!(ic, ::Val{2}, src1, src2, dest) = (ic[dest] = ic[src1] * ic[src2]; -1)
run_opcode!(ic, ::Val{3}, arg) = (ic[arg] = takeinput!(ic); -1)
run_opcode!(ic, ::Val{4}, arg) = (put!(ic.output, ic[arg]); -1)
run_opcode!(ic, ::Val{5}, pred, inst) = ic[pred] != 0 ? ic[inst] : -1
run_opcode!(ic, ::Val{6}, pred, inst) = ic[pred] == 0 ? ic[inst] : -1
run_opcode!(ic, ::Val{7}, arg1, arg2, dest) = (ic[dest] = Int(ic[arg1] < ic[arg2]); -1)
run_opcode!(ic, ::Val{8}, arg1, arg2, dest) = (ic[dest] = Int(ic[arg1] == ic[arg2]); -1)
run_opcode!(ic, ::Val{9}, arg) = (ic.relativeindex[] += ic[arg]; -1)

export compute!, IntComputer, getoutput, addinput!

end #OpCodes


text:
import PhysicsEngine

@testset "Grid" begin
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
"""
    postwalk(f, expr)
Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *after* they have been transformed by the walk. See also
`prewalk`.
"""
postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

"""
    prewalk(f, expr)
Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *before* they have been transformed by the walk, and the
walk will be applied to whatever `f` returns.
This makes `prewalk` somewhat prone to infinite loops; you probably want to try
`postwalk` first.
"""
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)

replace(ex, s, s′) = prewalk(x -> x == s ? s′ : x, ex)

"""
    inexpr(expr, x)
Simple expression match; will return `true` if the expression `x` can be found
inside `expr`.
    inexpr(:(2+2), 2) == true
"""
function inexpr(ex, x)
  result = false
  postwalk(ex) do y
    if y == x
      result = true
    end
    return y
  end
  return result
end