text:
module Sens
using LightXML, Printf
using DifferentialEquations, Statistics, GlobalSensitivity
using QuasiMonteCarlo
include("Constants.jl")
include("Utils.jl")
include("Batch.jl")
include("Plug.jl")
include("Cstr.jl")
include("Reactions.jl")
include("SurfaceReactions.jl")
struct  gsa_parameters
    monitor::Array{String,1}
    lower_pc::Float64
    upper_pc::Float64
    N::Int64
    gsa_srxn_ids::Array{Int64,1}
    gsa_srxn_constraint_ids::Array{Int64,1}
end
struct Mechanisms
    smd::MechanismDefinition
end

function gsa_xml(input_file::AbstractString, lib_dir::AbstractString)
    xmldoc = parse_file(input_file)
    xmlroot = root(xmldoc)

    # list of gasphase species
    gsa_species = get_collection_from_xml(xmlroot,"gsa_species")
    # lower bound for parameters (in percentages)
    lower_pc = get_value_from_xml(xmlroot,"gsa_lb")
    # upper bound for parameters (in percentages)
    upper_pc = get_value_from_xml(xmlroot,"gsa_ub")

    # Size of the sample
    sample_N = Int(get_value_from_xml(xmlroot,"gsa_N"))

    # get the reactions id which forms the parameter
    gsa_srxn_ids,gsa_srxn_constraint_ids = get_rxn_ids(xmlroot,"p_smech")

    gsa_p = gsa_parameters(gsa_species,lower_pc,upper_pc,sample_N,gsa_srxn_ids,gsa_srxn_constraint_ids)

    # get the reactor model
    model = get_text_from_xml(xmlroot,"gsa_model")

    model_root = root(parse_file(model))
    gsa_model = name(model_root)

    if lowercase(strip(gsa_model)) == "batch"
        params, prob, t_span = Batch.batch(model,lib_dir,true)
        gsa_problem(params, prob, t_span, gsa_p)
    elseif lowercase(strip(gsa_model)) == "plug"
        params, prob, t_span = Plug.plug(model,lib_dir, true)
        gsa_problem(params, prob, t_span, gsa_p)
    elseif lowercase(strip(gsa_model)) == "cstr"
        params, prob, t_span = Cstr.cstr(model,lib_dir,true)
        gsa_problem(params, prob, t_span, gsa_p)
    else
        println("Models yet to be implemented\n")
    end
end

"""
A common function for all reactor models to define the sensitivity problem
#   Usage
    gsa_problem(params, prob, gsa_p)
-   params : a tuple of parameters returned by the respective reactor model invoked
-   prob : the problem definition for the reactor model
-   gsa_p : struct of the type gsa_parameters
"""
function gsa_problem(params, prob, t_span, gsa_p::gsa_parameters)
    gsa_smech_params = Array{Float64,1}()
    state, thermo_obj, md, cp, o_stream = params
    n_surface_species = length(md.sm.species) # number of surface species
    n_gaspahse_species = length(thermo_obj.thermo_all) # number of gasphase species
    n_species = n_gaspahse_species + n_surface_species # total number of species
    # get the parameters from the surface reaction mechaism
    gsa_smech_parameter_map!(gsa_smech_params, gsa_p.gsa_srxn_ids, gsa_p.gsa_srxn_constraint_ids, md)
    println("Parameter size: ", length(gsa_smech_params))
    # create the lower limit for the parameters
    lb = gsa_smech_params .* (1-0.01*gsa_p.lower_pc)
    # create the upper limit for the parameters
    ub = gsa_smech_params .* (1+0.01*gsa_p.upper_pc)
    parameter_ratio = Array{Float64,1}()  # Array for storing the original parameter ratio (forward to reverse)
    if count(x->x>0, gsa_p.gsa_srxn_constraint_ids) > 0
        # get the initial parameter ratios if the reverse reaction constraints are specified
        initial_parameter_ratio!(parameter_ratio, gsa_p.gsa_srxn_ids, gsa_p.gsa_srxn_constraint_ids,md)
    end

    perform_gsa(params,prob, t_span, gsa_p, lb, ub, parameter_ratio)
end

function perform_gsa(params, prob, t_span, gsa_p::gsa_parameters,lb::Array{Float64,1}, ub::Array{Float64,1}, parameter_ratio::Array{Float64,1})
    state, thermo_obj, md, cp, o_stream = params
    time_points = collect(range(t_span[1],stop=t_span[2],length=50))
    #Get the name of all gasphase species
    species_list = Array{String,1}()
    for sp_thermo in thermo_obj.thermo_all
        push!(species_list,sp_thermo.name)
    end
    append!(species_list,md.sm.species)
    #get the id of species to be monitored from the gsa_species list
    monitor_species_ids = map(x->get_index(x,species_list),gsa_p.monitor)
    gsa_out = zeros(length(monitor_species_ids))
    func_eval_count = 0
    sens = function (gsa_params::Array{Float64})
        func_eval_count += 1
        println("Sample @ ", func_eval_count)
        SurfaceReactions.update_params!(md,gsa_params,gsa_p.gsa_srxn_ids,gsa_p.gsa_srxn_constraint_ids,parameter_ratio)
        params_updated = (state, thermo_obj, md, cp, o_stream)
        s_prob = remake(prob,p=params_updated)
        sol = solve(s_prob,alg_hints=[:stiff] , reltol=1e-6, abstol = 1e-8, saveat=time_points)
        for i in eachindex(monitor_species_ids)
            gsa_out[i] = mean(sol[monitor_species_ids[i],:])
        end
        gsa_out
    end
    sampler = SobolSample()
    A,B = QuasiMonteCarlo.generate_design_matrices(gsa_p.N,lb,ub,sampler)
    sobol_indices = gsa(sens,Sobol(order=[0,1,2]),A,B)

    # output file
    sobol_total = open("sobol_total.dat","w")
    sobol_first = open("sobol_first.dat","w")
    sobol_second = open("sobol_second.dat","w")
    header_data = [ "p"*string(k) for k in gsa_p.gsa_srxn_ids]
    create_header(sobol_total,["species"],header_data)
    create_header(sobol_first,["species"],header_data)
    create_header(sobol_second,["species"],header_data)

    for i in eachindex(monitor_species_ids)
        write_to_file(sobol_total,gsa_p.monitor[i],sobol_indices.ST[i,:])
        write_to_file(sobol_first,gsa_p.monitor[i],sobol_indices.ST[i,:])
        write_to_file(sobol_second,gsa_p.monitor[i],sobol_indices.ST[i,:])
    end

    close(sobol_total)
    close(sobol_first)
    close(sobol_second)

end

"""
A function to extract the parameters for sensitivity analysis from the xml input
    The function returnds two integer arrays. forward_rxn_ids contains the ids of the
    reactions whoes parameters will be considered as sensitivity parameters.
    reverse_rxn_ids contains the ids of the reactions aginst with the forward reaction
    parameters are correlated.
"""
function get_rxn_ids(xmlroot::XMLElement, xmltag::AbstractString)
    N = 10000
    p_mech = get_elements_by_tagname(xmlroot,xmltag)
    p_mech_content = ""
    if is_elementnode(p_mech[1])
        p_mech_content = content(p_mech[1])
    end
    p_mech_content_split = split(p_mech_content,",")
    forward_rxn_ids = Array{Int64,1}()
    reverse_rxn_ids = Array{Int64,1}()
    if strip(p_mech_content_split[1]) == ":"
        return [i for i in 1:N], [0 for i in 1:N]
    end
    for items in p_mech_content_split
        if occursin("=",items)
            rxn_sets = split(items,"=")
            forward_sets = rxn_sets[1]
            reverse_sets = rxn_sets[2]
            forward_sets = map(x->parse(Int64,x),split(forward_sets,":"))
            reverse_sets = map(x->parse(Int64,x),split(reverse_sets,":"))
            if length(forward_sets) > 1
                for ids in forward_sets[1]:forward_sets[2]
                    push!(forward_rxn_ids,ids)
                end
                for ids in reverse_sets[1]:reverse_sets[2]
                    push!(reverse_rxn_ids,ids)
                end
            else
                push!(forward_rxn_ids,forward_sets[1])
                push!(reverse_rxn_ids,reverse_sets[1])
            end
        else
            forward_sets = map(x->parse(Int64,x),split(items,":"))
            if length(forward_sets) > 1
                for ids in forward_sets[1]:forward_sets[2]
                    push!(forward_rxn_ids,ids)
                    push!(reverse_rxn_ids,0)
                end
            else
                push!(forward_rxn_ids,forward_sets[1])
                push!(reverse_rxn_ids,0)
            end
        end
    end
    return (forward_rxn_ids, reverse_rxn_ids)
end

"""
This function returns the sticking coefficient or the pre-exponent factor
    corresponding to the reaction ids present in gsa_rxn_ids
#   Usage
    gsa_smech_parameter_map!(gsa_smech_params, gsa_srxn_ids, gsa_srxn_constraint_ids, md )
-   gsa_smech_params : vector containing the parameters extracted from the inoput mechanism (output)
-   gsa_srxn_ids : reaction ids whoes parameters needs to be extracted from the mechanism (input)
-   gsa_srxn_constraint_ids : id of the reactions which needs to be considered for reverse reaction
-   md: MechanismDefinition
"""
function gsa_smech_parameter_map!(gsa_smech_params::Array{Float64,1}, gsa_srxn_ids::Array{Int64,1}, gsa_srxn_constraint_ids::Array{Int64,1}, md )

    function populate(rxn)
        if Int(rxn.type) == Int(stick)
            push!(gsa_smech_params,rxn.params.s)
        else
            push!(gsa_smech_params,rxn.params.k0)
        end
    end
    # if the length of gsa_srxn_ids is very high then it means all reactions parameters need to be considered
    if length(gsa_srxn_ids) >= length(md.sm.reactions)
        for rxn in md.sm.reactions
            populate(rxn)
        end
        gsa_srxn_ids_copy = copy(gsa_srxn_ids)
        empty!(gsa_srxn_ids)
        gsa_srxn_constraint_ids_copy = copy(gsa_srxn_constraint_ids)
        empty!(gsa_srxn_constraint_ids)
        for i in 1:length(md.sm.reactions)
            push!(gsa_srxn_ids,gsa_srxn_ids_copy[i])
            push!(gsa_srxn_constraint_ids,gsa_srxn_constraint_ids_copy[i])
        end
        empty!(gsa_srxn_ids_copy)
        empty!(gsa_srxn_constraint_ids_copy)
    else
        for rxn in md.sm.reactions
            if rxn.id in gsa_srxn_ids
                populate(rxn)
            end
        end
    end

    return gsa_smech_params
end

"""
A function to calculate the intial ratio of pre-exponential factors
"""
function initial_parameter_ratio!(p_ratio::Array{Float64,1},rxn_ids::Array{Int64,1},constraints_ids::Array{Int64,1}, md)
    for i in eachindex(rxn_ids)
        if constraints_ids[i] != 0
            if Int(md.sm.reactions[rxn_ids[i]].type) == Int(stick)
                ratio = md.sm.reactions[rxn_ids[i]].params.k0/md.sm.reactions[constraints_ids[i]].params.k0
                push!(p_ratio,ratio)
            else
                ratio = md.sm.reactions[rxn_ids[i]].params.k0/md.sm.reactions[constraints_ids[i]].params.k0
                push!(p_ratio,ratio)
            end
        else
            push!(p_ratio,0.0)
        end
    end
end

#end of module
end


text:
module sNMF_VB

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions

export standard_VB, online_VB

function standard_VB(X::Array{ℜ,2}, R::Ƶ ; μ::ℜ=1.0, γ::ℜ=0.1, EPOCHS::Ƶ=1, ϵ::ℜ=1e-16) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ = size(X,1)

    a::ℜ, b::ℜ = I*I*γ, γ/μ
    α_R, α_IR = fill(a/R,R), fill(2.0*a/(I*R),I,R)
    S₊ = sum(X)

    log_λ = digamma(S₊ + a) - log(b + 1)
    log_θ_R = log.(rand(Dirichlet(α_R.+0.1)))
    log_θ_IR = reshape([log(θ_ir) for r=1:R for θ_ir=rand(Dirichlet(α_IR[:,r].+0.1))],I,R)

    S_R, S_IR, S⁺::ℜ = zeros(ℜ,R), zeros(ℜ,I,R), 0.0
    s, p, log_ρ = Array{ℜ}(undef,R), Array{ℜ}(undef,R), Array{ℜ}(undef,R)

    ELBO::Array{ℜ} = zeros(ℜ,EPOCHS)

    X_full = similar(X)
    for eph=1:EPOCHS
        S_R .= 0.0
        S_IR .= 0.0
        S⁺ = 0.0

        for i₂=1:I, i₁=1:I #order of traversal is important
            log_ρ .= log_λ .+ log_θ_R .+ log_θ_IR[i₁,:] .+ log_θ_IR[i₂,:]
            log_ρ₊ = logsumexp(log_ρ)
            p .= exp.(log_ρ .- log_ρ₊)

            X_full[i₁,i₂] = isnan(X[i₁,i₂]) ? exp(log_ρ₊) : X[i₁,i₂]
            s .= X_full[i₁,i₂] .* p

            S_IR[i₁,:] .+= s
            S_IR[i₂,:] .+= s
            S_R .+= s
            S⁺ += X_full[i₁,i₂]

            ELBO[eph] += X[i₁,i₂]*log_ρ₊ - lgamma(X[i₁,i₂] + 1.0)
        end

        ELBO[eph] += a*log(b) - (a + S⁺)*log(b+1.0) - sum(lgamma.(α_R)) + sum(lgamma.(α_R .+ S_R))
        ELBO[eph] += sum(lgamma.(2.0 .* α_R)) - sum(lgamma.(α_IR))
        ELBO[eph] += sum(lgamma.(α_IR .+ S_IR)) - sum(lgamma.(2.0 .* (α_R .+ S_R)))
        ELBO[eph] -= S⁺ * log_λ + sum(S_R .* log_θ_R) + sum(S_IR .* log_θ_IR)

        log_λ = digamma(S⁺+a) - log(b + 1.0)
        log_θ_R .= digamma.(S_R .+ α_R) .- digamma(S⁺+a)
        log_θ_IR .= digamma.(S_IR .+ α_IR) .- digamma.(S_R.+α_R)'

    end
    return ELBO, X_full, log_λ, (log_θ_R, log_θ_IR)
end

end

ext: jl

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
using DataFrames, DataFramesMeta, CSV, CodecZlib
using Plots, StatsPlots, Measures
using BioSequences, FASTX, GenomicFeatures, GenomicIntersections, GenomeFragments
using MotifScanner
using ProgressMeter
using Statistics, StatsBase

## Some packages are loaded within included files below
include("mutations.jl")
include("plots.jl")
include("motifs.jl")
include("regions.jl")
include("genomics.jl")
include("conformation.jl")
include("expression.jl")

"""
    loaddata(loadalignments=true, loadhic=true; coolerfile="")

    Loads data necessary for analysis and plotting. This repo contains all source data necessary to regenerate the figures calculated from original alignment files and cooler file.
    Should you wish to recalculate source files, then include alignment files in data/chip_atac_datasets.tsv (.bin format see https://github.com/owensnick/GenomeFragments.jl) and
    the cooler file supplied as an argument to this function. Then set loadalignments=true and loadhic=true.

"""
function loaddata(loadalignments=false, loadhic=false; coolerfile="")

    ### mutations and region
    mutations = loadmutations()
    region = loadregregion()
    annotateregion!(mutations, region)
    distinctmutations, refseq, altseqs = distinctmuts(mutations, region);
    hk1coords = hk1_locus_coords()
    gtf = loadgtf()
    genecoords = transcriptcoords(gtf);

    ### chip/atac datasets
    datasets = loaddatasets();

    datagroup = combine(groupby(datasets, [:Study, :DataType, :Cell, :Target]), nrow => :count, :Index => Ref => :Inds);
    datagroup.Label = @with datagroup string.(:Study, "_", :DataType, "_", :Cell, "_", :Target);
    sort!(datagroup, [:Study, :Target, :Cell]);

    if loadalignments
        if any(ismissing, datasets.SignalFile)
            @warn "Signal files missing from data/chip_atac_datasets.tsv. Proceeding without loading."
            FM = nothing
        else
            FM = load_frag_matrix.(datasets.SignalFile);
        end
    else
        FM = nothing
    end

    ### hicdata
    hicdata = loadhicdata();

    ### Expression datasets
    gtex = loadgtexgct();
    betadiffstage = load_beta_diff_geusz_xie();
    betadiffsc = load_betadiff_sc();
    islet_scrnaseq = loadisletscrna();
    betadifflate = load_betadiff_scrnaseq();

    ### cooler
    if loadhic

        if !isfile(replace(coolerfile, r"::resolutions\/[0-9]*$" => "")) ## strip off cooler resolution
            @warn "Cooler file $coolerfile not found. Proceeding without loading."
            isletcooler = nothing
        else
            isletcooler = loadcoolerfile(coolerfile)
        end

    else
        isletcooler = nothing
    end
    hicdata = quantisletloops(hicdata, isletcooler) ## quantify HIC data in islet loops

    mutations, region, distinctmutations, refseq, altseqs, hk1coords, genecoords, datasets, datagroup, FM, hicdata, gtex, betadiffstage, betadiffsc, islet_scrnaseq, betadifflate, isletcooler
end

#### Displaying wide tables, primary use in jupyter notebooks
function showwide(table)
    c = ENV["COLUMNS"]
    ENV["COLUMNS"] = "10000"
    display(table)
    ENV["COLUMNS"] = c;
    nothing;
end

#### Displaying wide tables and long tables, returns a function that will do displaying
function showwl(lines=200, nc=10000)

    table -> begin
        l = ENV["LINES"]
        c = ENV["COLUMNS"]
        ENV["LINES"] = string(lines)
        ENV["COLUMNS"] = string(nc)
        display(table)
        ENV["LINES"] = l;
        ENV["COLUMNS"] = c;
        nothing;
    end
end

### function to get base dir of project whether working in src or notebooks
function getprojectdir()
    d = pwd()
    if (basename(d) == "notebooks") || (basename(d) == "src")
        return dirname(d)
    else
        return d ## assume working in project dir
    end
end


text:
filename = [
    "merge_sort",
    "quick_sort",
    "binary_search",
    "tree_traverse",
    "binary_search_tree",
    "bfs",
    "dfs",
    "topological_sort",
    "monte_carlo_integral",
    "mcm",
    "mst"
]

for f in filename
    teststr = "# Test '$(f)' #"
    len = length(teststr)
    teststr = "\n$("#"^len)\n$(teststr)\n$("#"^len)\n"
    println()
    @info teststr

    include("$(f).jl")
end


text:
export get_tempdir,
       get_tempfilepath,
       util_open,
       util_close,
       load_data_file,
       dict_remap,
       NativeEncodingToUnicode,
       PDFEncodingToUnicode,
       utilPrintOpenFiles,
       to_uint8,
       to_uint16

using Pkg

# Non-exported variable that retains the current tempdir location.
# This shoould be only accessed by get_tempdir() only

CURRENT_TMPDIR=""

"""
On every launch of the application a new directory needs to be created
that will have the path of the new directory.
"""
function get_tempdir()
    global CURRENT_TMPDIR
    if (CURRENT_TMPDIR == "" )||!isdir(CURRENT_TMPDIR)
        CURRENT_TMPDIR= abspath(mktempdir())
    end
    return CURRENT_TMPDIR
end

get_tempfilepath()=mktemp(get_tempdir())
util_open(filename, mode)=open(filename, mode)
util_close(handle)=close(handle)
utilPrintOpenFiles()=[]

import Base: zero
zero(::Type{Char}) = Char(0x00)

using DelimitedFiles

function load_data_file(filename)
    path = joinpath(@__DIR__, "..", "data", filename)
    return readdlm(path, ',', String, '\n',  comments=true, comment_char='#')
end

function dict_remap(ab, bc)
    d = Dict()
    for (a, b) in ab
        c = get(bc, b, zero(valtype(bc)))
        d[a] = c
    end
    return d
end

@static if VERSION >= v"0.7-"
    to_uint16(x) = parse(UInt16, x, base=16)
    to_uint8(x)  = parse(UInt8,  x, base=8)
else
    to_uint16(x) = parse(UInt16, x, 16)
    to_uint8(x)  = parse(UInt8,  x, 8)
end

const PDFEncoding_to_Unicode = begin
    d = Dict()
    m = load_data_file("pdf-doc-encoding.txt")
    map(m[:,3], m[:,4]) do x, y
        e = to_uint8(x)
        u = (y != "") ? Char(to_uint16(y)) : Char(e)
        d[e] = u
    end
    d
end

function NativeEncodingToUnicode(barr::Vector{UInt8}, mapping::Dict)
    l = length(barr)
    carr = similar(barr, Char)
    for i = 1:l
        carr[i] = get(mapping, barr[i], zero(Char))
    end
    return carr
end

PDFEncodingToUnicode(barr::Vector{UInt8}) =
    NativeEncodingToUnicode(barr, PDFEncoding_to_Unicode)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module SortPerf
import Perftests: @perf, meta

using Base.Sort
import Base.Sort: QuickSort, MergeSort, InsertionSort

Pkg.add("SortingAlgorithms")
using SortingAlgorithms #Provides the other sorting algorithms

sorts = [InsertionSort, QuickSort, MergeSort, HeapSort, RadixSort, TimSort]

randstr_fn!(str_len::Int) = d -> (for i = 1:length(d); d[i] = randstring(str_len); end; d)
randint_fn!(m::Int) = d -> rand!(d, 1:m)

# Run the full gamut of tests
for (T, typename, randfn!) in [(Int, string(Int), randint_fn!(10)),
                                (Float64, string(Float64), rand!),
                                (AbstractString, "String_05", randstr_fn!(5)),
                                (AbstractString, "String_10", randstr_fn!(10))]
    for logsize = 6:2:18
        size = 2^logsize
        for s in sorts
            if s == RadixSort && T == AbstractString continue end      #Radix sort not implemented
            if s == InsertionSort && logsize >=14 continue end #Too slow

            data = Array(T, size)
            gc()

            # Note; we're using an anonymous function here so that
            perfsort = (data) -> sort(data, alg=s)

            # Randomize the data, measure the sorting, then measure the sorting again
            typesize = "$(typename)_$(logsize)"
            sname = string(s)[1:end-5]
            sname = sname[rsearch(sname,'.')+1:end]

            randfn!(data)
            @perf perfsort(data) meta(sname, "$(typesize)_random", "Sorting of randomized data")
            sort!(data)
            @perf perfsort(data) meta(sname, "$(typesize)_sorted", "Sorting of sorted data")

            # Measure it after reversing the data
            rdata = reverse(data)
            @perf perfsort(rdata) meta(sname, "$(typesize)_reverse", "Sorting of reverse-sorted data")

            ## Sorted with 3 exchanges
            for i = 1:3
                n1 = rand(1:size)
                n2 = rand(1:size)
                data[n1], data[n2] = data[n2], data[n1]
            end
            @perf perfsort(data) meta(sname, "$(typesize)_3exchanges", "Sorting of sorted data with three exchanges")

            ## Sorted with 10 unsorted values appended
            sort!(data)
            data[end-9:end] = randfn!(Array(T,10))
            @perf perfsort(data) meta(sname, "$(typesize)_appended", "Sorting of data with 10 appended unsorted values")

            ## Random data with 4 unique values
            data4 = data[rand(1:4,size)]
            @perf perfsort(data4) meta(sname, "$(typesize)_4unique", "Sorting of data with only 4 unique values")

            ## All values equal
            data1 = data[ones(Int, size)]
            @perf perfsort(data1) meta(sname, "$(typesize)_allequal", "Sorting of completely homogenous data")

            ## QuickSort median killer
            if s == QuickSort && logsize > 16; continue; end  # too slow!

            data = data[1:size>>1]
            data = sort!(data)
            data = vcat(reverse(data), data)
            @perf perfsort(data) meta(sname, "$(typesize)_qsortkiller", "Sorting of qsort worst case")
        end
    end
end

end # module


text:
count = 0

function sleep_and_print(n)
    global count
    sleep(n)
    count += 1
    println(n)
end

a = [5, 3, 6, 2, 3, 2, 2];
for i=1:length(a)
    @async sleep_and_print(a[i])
end

# wait for all threads to complete
while count != length(a)
    sleep(1)
end

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
# DESCRIPTION
# Simple model of a helicopter using static aerodynamic model for rotor,
# tail rotor, horzontal and vertical tail, and fuselage.
#
# INPUT
# - state: state vector
# - ctrls: stick control inputs
#
# OUTPUT
# - dstate: state derivative
# - FM: total forces and moments
#
# ------------------------------------------------------------------------------

function SimpleHel!(state,sctrls)

      # unpack states
      u=state[1]
      v=state[2]
      w=state[3]
      p=state[4]
      q=state[5]
      r=state[6]
      phi=state[7]
      theta=state[8]
      psi=state[9]
      x=state[10]
      y=state[11]
      z=state[12]
      # convert to stick inputs to swashplate and TR collective inputs
      ctrls = control_mixing(sctrls)
      # unpack controls
      theta1c=ctrls[1]*pi/180
      theta1s=ctrls[2]*pi/180
      theta0=ctrls[3]*pi/180
      theta0tr=ctrls[4]*pi/180

      # ----------------------------- MAIN ROTOR -------------------------------

      # hub velocities
      Vh=[u;v;w]+cross([p;q;r],rMR)
      uh=Vh[1]
      vh=Vh[2]
      wh=Vh[3]
      # advance ratios
      Vinplane=sqrt(uh^2+vh^2)
      mu=Vinplane/(OMEGA*RADIUS)
      muz=wh/(OMEGA*RADIUS)
      mu2=mu*mu
      # get wind axis transformation
      psiw=atan2(vh,uh)
      cpsiw=cos(psiw)
      spsiw=sin(psiw)
      Th2w=[ cpsiw spsiw
            -spsiw cpsiw]
      Tw2h=Th2w'
      # non-dimensional angular rates in wind axes
      temp=Th2w*[p; q]
      phw=temp[1]/OMEGA
      qhw=temp[2]/OMEGA
      # cyclic pitch translated to wind axes
      temp=Th2w*[theta1s;theta1c]
      theta1sw=temp[1]
      theta1cw=temp[2]
      # CT and lambda0 iteration, solves for inflow and thrust coefficient
      # first guess
      lambda0=0.05
      # set tolerance greater than error
      delta_lambda=1
      iter=0
      function inflow_MR()
            while (abs(delta_lambda)>1e-10 && iter<200)
                  CT=0.5*A0*SOL*((1/3+0.5*mu2)*theta0+0.5*mu*(theta1sw+0.5*phw)+
                     0.5*(muz-lambda0)+0.25*(1+mu2)*TWIST)
                  lamtot2=mu2+(lambda0-muz)^2
                  delta_lambda=-(2.0*lambda0*sqrt(lamtot2)-CT)*lamtot2/(2.0*lamtot2^1.5+0.25*
                               A0*SOL*lamtot2-CT*(muz-lambda0))
                  lambda0=lambda0+0.5*delta_lambda
            iter=iter+1
            end
            return CT
      end
      CT=inflow_MR()
      # thrust
      T=CT*RHO*OMEGA^2*pi*RADIUS^4
      # uniform inflow assumption in this model
      lambda1cw=0.0
      lambda1sw=0.0
      # quasi-Steady flapping in wind axes (Padfield p. 107)
      Abt = [-8/3*mu*(1+0.5*mu2)    -2*mu*(1+0.5*mu2)        -(1+2*mu^2)      0.0
             -GAMMA/6*mu*(1+0.5*mu2) -2*GAMMA*mu/15*(1+mu2/3) -2*GAMMA*mu2/9.0 (1-0.5*mu2*mu2)]
      Abl = [-2*mu*(1+0.5*mu2)          (1+0.5*mu^2)                    0.0
             -2*GAMMA*mu/9.0*(1-0.5*mu2)  GAMMA/9.0*mu^2+0.5*GAMMA/9.0*mu2  -(1.0-0.5*mu2) ]
      Abo = [-(1+0.5*mu^2)                          16/GAMMA*(1+0.5*mu2)
             16/GAMMA*(1.0-0.5*mu2)+GAMMA/9.0*mu^2  (1.0-0.5*mu2)]
      #
      temp=Abt*[theta0; TWIST; theta1sw; theta1cw]+Abl*
               [muz-lambda0; lambda1sw; lambda1cw]+Abo*[phw; qhw]
      beta1cw=temp[1]
      beta1sw=temp[2]
      # convert back to hub system
      temp=Tw2h*[beta1sw; beta1cw]
      beta1s=temp[1]
      beta1c=temp[2]
      # torque
      delta=DELTA0+DELTA2*CT*CT
      CQ=-(muz-lambda0)*CT+0.125*delta*SOL*(1+(7/3)*mu^2)
      Q=CQ*RHO*OMEGA^2*pi*RADIUS^5
      # translate rotor forces and moments to CG
      Xr=T*beta1c
      Yr=-T*beta1s
      Zr=-T
      MvecMR=[0.0; 0.0; Q]+cross(rMR, [Xr; Yr; Zr])
      Lr=MvecMR[1]
      Mr=MvecMR[2]
      Nr=MvecMR[3]

      # ------------------------------ FUSELAGE --------------------------------

      # velocity vector
      Vf=[u; v; w]
      # drag force vector
      Ff=-0.5*RHO*sqrt(Vf'*Vf)*Vf*CDSF
      # transform to body forces
      Xf=Ff[1]
      Yf=Ff[2]
      Zf=Ff[3]
      MvecF=cross(rFS,[Xf;Yf;Zf])
      Lf=MvecF[1]
      Mf=MvecF[2]
      Nf=MvecF[3]

      # ----------------------- HORIZONTAL STABILIZER --------------------------

      # local angle of attack
      Vht=[u; v; w]+cross([p; q; r], rHT)
      uht=Vht[1]
      wht=Vht[3]
      alphaht=atan2(wht,uht)
      # lift force
      CLht=AHT*alphaht
      LIFTht=0.5*RHO*(Vht'*Vht)*SHT*CLht
      # transform to body forces at CG
      Xht=0.0
      Yht=0.0
      Zht=-LIFTht
      MvecHT=cross(rHT,[Xht;Yht;Zht])
      Lht=MvecHT[1]
      Mht=MvecHT[2]
      Nht=MvecHT[3]

      # ----------------------- VERTICAL STABILIZER --------------------------

      # local angle of attack
      Vvt=[u; v; w]+cross([p; q; r], rVT)
      uvt=Vvt[1]
      vvt=Vvt[2]
      alphavt=atan2(vvt,uvt)
      # lift force
      CLvt=AVT*alphavt
      LIFTvt=0.5*RHO*(Vvt'*Vvt)*SVT*CLvt
      # transform to body forces at CG
      Xvt=0.0
      Yvt=LIFTvt
      Zvt=0.0
      MvecVT=cross(rVT,[Xvt;Yvt;Zvt])
      Lvt=MvecVT[1]
      Mvt=MvecVT[2]
      Nvt=MvecVT[3]

      # ----------------------------- TAIL ROTOR -------------------------------

      # local velocities at tail rotor hub in body CS
      Vtr=[u;v;w]+cross([p;q;r],rTR)
      # transfrom to TR coordinate system
      Ttrg=[1.0  0.0 0.0
            0.0  0.0 1.0
            0.0 -1.0 0.0]
      Vtr_tr=Ttrg*Vtr
      utr=Vtr_tr[1]
      vtr=Vtr_tr[2]
      wtr=Vtr_tr[3]
      mutr=sqrt(utr^2+vtr^2)/(OMEGTR*RTR)
      muztr=wtr/(OMEGTR*RTR)
      # CT and lambda0 iteration for tail rotor
      # first guess
      lambda0tr=0.05
      #set tolerance greater than error
      delta_lambda=1
      iter=0
      function inflow_TR()
            while(abs(delta_lambda)>1e-10 && iter<200)
                  CTtr=0.5*A0TR*SOLTR*((1/3+0.5*mutr^2)*theta0tr+0.5*(muztr-lambda0tr)+
                       0.25*(1+mutr^2)*TWISTTR)
                  lamtot2=mutr^2+(lambda0tr-muztr)^2
                  delta_lambda=-(2.0*lambda0tr*sqrt(lamtot2)-CTtr)*lamtot2/(2.0*lamtot2^1.5+
                               0.25*A0TR*SOLTR*lamtot2-CTtr*(muztr-lambda0tr))
                  lambda0tr=lambda0tr+0.5*delta_lambda
                  iter=iter+1
            end
            return CTtr
      end
      CTtr=inflow_TR()
      # tail rotor thrust
      Ttr=CTtr*RHO*OMEGTR^2*pi*RTR^4
      # force vector in body coordinates
      Ftr=Ttrg'*[0.0;0.0;-Ttr]
      Xtr=Ftr[1]
      Ytr=Ftr[2]
      Ztr=Ftr[3]
      # resolve force and moment about CG
      MvecTR=cross(rTR,Ftr)
      Ltr=MvecTR[1]
      Mtr=MvecTR[2]
      Ntr=MvecTR[3]

      # ------------------------ EQUATIONS OF MOTION ---------------------------

      # sum forces and moments
      X=Xr+Xf+Xtr+Xht+Xvt
      Y=Yr+Yf+Ytr+Yht+Yvt
      Z=Zr+Zf+Ztr+Zht+Zvt
      L=Lr+Lf+Ltr+Lht+Lvt
      M=Mr+Mf+Mtr+Mht+Mvt
      N=Nr+Nf+Ntr+Nht+Nvt
      # total forces and moments
      FM=zeros(6,1)
      FM=[X;Y;Z;L;M;N]

      # EQUATIONS OF MOTION
      # trigonometric functions
      cphi=cos(phi)
      sphi=sin(phi)
      cthe=cos(theta)
      sthe=sin(theta)
      cpsi=cos(psi)
      spsi=sin(psi)
      # declare state derivative vector
      dstate=zeros(12,1)
      # equations of motion
      # translational dynamics
      dstate[1]=X/mass-g*sthe-q*w+r*v
      dstate[2]=Y/mass+g*cthe*sphi-r*u+p*w
      dstate[3]=Z/mass+g*cthe*cphi-p*v+q*u
      # rotational dynamics
      gam=Ix*Iz-Ixz^2
      dstate[4]=(Iz*L+Ixz*N+Ixz*(Ix-Iy+Iz)*p*q-(Iz^2-Iy*Iz+Ixz^2)*q*r)/gam
      dstate[5]=(M+(Iz-Ix)*p*r-Ixz*(p^2-r^2))/Iy
      dstate[6]=(Ix*N+Ixz*L-Ixz*(Ix-Iy+Iz)*q*r+(Ix^2-Ix*Iy+Ixz^2)*p*q)/gam
      # rotational kinematics
      dstate[7]=p+q*sphi*sthe/cthe+r*cphi*sthe/cthe
      dstate[8]=q*cphi-r*sphi
      dstate[9]=q*sphi/cthe+r*cphi/cthe
      # position
      dstate[10]=u*cthe*cpsi+v*(sphi*sthe*cpsi-cphi*spsi)+w*(cphi*sthe*cpsi+sphi*spsi)
      dstate[11]=u*cthe*spsi+v*(sphi*sthe*spsi+cphi*cpsi)+w*(cphi*sthe*spsi-sphi*cpsi)
      dstate[12]=-u*sthe+v*sphi*cthe+w*cphi*cthe
      #
      return dstate, FM
end

max_stars_repo_name: umbertosaetti/J-SimpleHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
# Differentiating tensor-product polynomials

function chebyshev_derivative(weights::Array{T,1},x::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain)

  poly_deriv_1 = chebyshev_polynomial_derivative(order[1],[x1])

  evaluated_derivative = zero(T)

  for i = 1:order[1]+1

    evaluated_derivative += (2.0/(domain[1]-domain[2]))*poly_deriv_1[i]*weights[i]

  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,2},x::Array{T,1},order::Array{S,1},domain=[ones(1,2);-ones(1,2)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])

  polynomial_1_data = zeros(2,order[1]+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order[1],x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order[1],x1)

  polynomial_2_data = zeros(2,order[2]+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order[2],x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order[2],x2)

  evaluated_derivative = zeros(1,2)

  for j = 1:order[2]+1
    for i = 1:order[1]+1

      evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*weights[i,j]
      evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*weights[i,j]

    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,3},x::Array{T,1},order::Array{S,1},domain=[ones(1,3);-ones(1,3)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])

  polynomial_1_data = zeros(2,order[1]+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order[1],x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order[1],x1)

  polynomial_2_data = zeros(2,order[2]+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order[2],x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order[2],x2)

  polynomial_3_data = zeros(2,order[3]+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order[3],x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order[3],x3)

  evaluated_derivative = zeros(1,3)

  for k = 1:order[3]+1
    for j = 1:order[2]+1
      for i = 1:order[1]+1

        evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*weights[i,j,k]
        evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[1,k]*weights[i,j,k]
        evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*weights[i,j,k]

      end
    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,4},x::Array{T,1},order::Array{S,1},domain=[ones(1,4);-ones(1,4)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])

  polynomial_1_data = zeros(2,order[1]+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order[1],x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order[1],x1)

  polynomial_2_data = zeros(2,order[2]+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order[2],x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order[2],x2)

  polynomial_3_data = zeros(2,order[3]+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order[3],x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order[3],x3)

  polynomial_4_data = zeros(2,order[4]+1)
  polynomial_4_data[1,:] = chebyshev_polynomial(order[4],x4)
  polynomial_4_data[2,:] = chebyshev_polynomial_derivative(order[4],x4)

  evaluated_derivative = zeros(1,4)

  for l = 1:order[4]+1
    for k = 1:order[3]+1
      for j = 1:order[2]+1
        for i = 1:order[1]+1

          evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*weights[i,j,k,l]
          evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[2,k]*polynomial_4_data[2,l]*weights[i,j,k,l]
          evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*polynomial_4_data[1,l]*weights[i,j,k,l]
          evaluated_derivative[4] += (2.0/(domain[1,4]-domain[2,4]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[2,l]*weights[i,j,k,l]

  			 end
			end
    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,5},x::Array{T,1},order::Array{S,1},domain=[ones(1,5);-ones(1,5)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])

  polynomial_1_data = zeros(2,order[1]+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order[1],x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order[1],x1)

  polynomial_2_data = zeros(2,order[2]+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order[2],x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order[2],x2)

  polynomial_3_data = zeros(2,order[3]+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order[3],x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order[3],x3)

  polynomial_4_data = zeros(2,order[4]+1)
  polynomial_4_data[1,:] = chebyshev_polynomial(order[4],x4)
  polynomial_4_data[2,:] = chebyshev_polynomial_derivative(order[4],x4)

  polynomial_5_data = zeros(2,order[5]+1)
  polynomial_5_data[1,:] = chebyshev_polynomial(order[5],x5)
  polynomial_5_data[2,:] = chebyshev_polynomial_derivative(order[5],x5)

  evaluated_derivative = zeros(1,5)

  for m = 1:order[5]+1
	  for l = 1:order[4]+1
      for k = 1:order[3]+1
        for j = 1:order[2]+1
          for i = 1:order[1]+1

            evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
            evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
            evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
            evaluated_derivative[4] += (2.0/(domain[1,4]-domain[2,4]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[2,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
            evaluated_derivative[5] += (2.0/(domain[1,5]-domain[2,5]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[2,m]*weights[i,j,k,l,m]

					end
				end
      end
    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,6},x::Array{T,1},order::Array{S,1},domain=[ones(1,6);-ones(1,6)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])
  x6 = normalize_node(x[6],domain[:,6])

  polynomial_1_data = zeros(2,order[1]+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order[1],x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order[1],x1)

  polynomial_2_data = zeros(2,order[2]+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order[2],x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order[2],x2)

  polynomial_3_data = zeros(2,order[3]+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order[3],x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order[3],x3)

  polynomial_4_data = zeros(2,order[4]+1)
  polynomial_4_data[1,:] = chebyshev_polynomial(order[4],x4)
  polynomial_4_data[2,:] = chebyshev_polynomial_derivative(order[4],x4)

  polynomial_5_data = zeros(2,order[5]+1)
  polynomial_5_data[1,:] = chebyshev_polynomial(order[5],x5)
  polynomial_5_data[2,:] = chebyshev_polynomial_derivative(order[5],x5)

  polynomial_6_data = zeros(2,order[6]+1)
  polynomial_6_data[1,:] = chebyshev_polynomial(order[6],x6)
  polynomial_6_data[2,:] = chebyshev_polynomial_derivative(order[6],x6)

  evaluated_derivative = zeros(1,6)

  for n = 1:order[6]+1
	  for m = 1:order[5]+1
	    for l = 1:order[4]+1
        for k = 1:order[3]+1
          for j = 1:order[2]+1
            for i = 1:order[1]+1

              evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
              evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
              evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
              evaluated_derivative[4] += (2.0/(domain[1,4]-domain[2,4]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[2,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
              evaluated_derivative[5] += (2.0/(domain[1,5]-domain[2,5]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[2,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
              evaluated_derivative[6] += (2.0/(domain[1,6]-domain[2,6]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[2,n]*weights[i,j,k,l,m,n]

            end
          end
				end
      end
    end
  end

  return evaluated_derivative

end

# Differentiating complete polynomials

function chebyshev_derivative(weights::Array{T,1},x::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain)

  poly_deriv_1 = chebyshev_polynomial_derivative(order,[x1])

  evaluated_derivative = zero(T)

  for i = 1:order+1

    evaluated_derivative += (2.0/(domain[1,1]-domain[2,1]))*poly_deriv_1[i]*weights[i]

  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,2},x::Array{T,1},order::S,domain=[ones(1,2);-ones(1,2)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])

  polynomial_1_data = zeros(2,order+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order,x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order,x1)

  polynomial_2_data = zeros(2,order+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order,x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order,x2)

  evaluated_derivative = zeros(1,2)

  for j = 1:order+1
    for i = 1:order+1

      if (i+j <= order+2)

        evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*weights[i,j]
        evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*weights[i,j]

      end

    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,3},x::Array{T,1},order::S,domain=[ones(1,3);-ones(1,3)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])

  polynomial_1_data = zeros(2,order+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order,x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order,x1)

  polynomial_2_data = zeros(2,order+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order,x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order,x2)

  polynomial_3_data = zeros(2,order+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order,x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order,x3)

  evaluated_derivative = zeros(1,3)

  for k = 1:order+1
    for j = 1:order+1
      for i = 1:order+1

        if (i+j+k <= order+3)

          evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*weights[i,j,k]
          evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[1,k]*weights[i,j,k]
          evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*weights[i,j,k]

        end

      end
    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,4},x::Array{T,1},order::S,domain=[ones(1,4);-ones(1,4)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])

  polynomial_1_data = zeros(2,order+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order,x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order,x1)

  polynomial_2_data = zeros(2,order+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order,x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order,x2)

  polynomial_3_data = zeros(2,order+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order,x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order,x3)

  polynomial_4_data = zeros(2,order+1)
  polynomial_4_data[1,:] = chebyshev_polynomial(order,x4)
  polynomial_4_data[2,:] = chebyshev_polynomial_derivative(order,x4)

  evaluated_derivative = zeros(1,4)

  for l = 1:order+1
    for k = 1:order+1
      for j = 1:order+1
        for i = 1:order+1

          if (i+j+k+l <= order+4)

            evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*weights[i,j,k,l]
            evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[2,k]*polynomial_4_data[2,l]*weights[i,j,k,l]
            evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*polynomial_4_data[1,l]*weights[i,j,k,l]
            evaluated_derivative[4] += (2.0/(domain[1,4]-domain[2,4]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[2,l]*weights[i,j,k,l]

          end

        end
			end
    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,5},x::Array{T,1},order::S,domain=[ones(1,5);-ones(1,5)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])

  polynomial_1_data = zeros(2,order+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order,x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order,x1)

  polynomial_2_data = zeros(2,order+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order,x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order,x2)

  polynomial_3_data = zeros(2,order+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order,x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order,x3)

  polynomial_4_data = zeros(2,order+1)
  polynomial_4_data[1,:] = chebyshev_polynomial(order,x4)
  polynomial_4_data[2,:] = chebyshev_polynomial_derivative(order,x4)

  polynomial_5_data = zeros(2,order+1)
  polynomial_5_data[1,:] = chebyshev_polynomial(order,x5)
  polynomial_5_data[2,:] = chebyshev_polynomial_derivative(order,x5)

  evaluated_derivative = zeros(1,5)

  for m = 1:order+1
	  for l = 1:order+1
      for k = 1:order+1
        for j = 1:order+1
          for i = 1:order+1

            if (i+j+k+l+m <= order+5)

              evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
              evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
              evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
              evaluated_derivative[4] += (2.0/(domain[1,4]-domain[2,4]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[2,l]*polynomial_5_data[1,m]*weights[i,j,k,l,m]
              evaluated_derivative[5] += (2.0/(domain[1,5]-domain[2,5]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[2,m]*weights[i,j,k,l,m]

            end

					end
				end
      end
    end
  end

  return evaluated_derivative

end

function chebyshev_derivative(weights::Array{T,6},x::Array{T,1},order::S,domain=[ones(1,6);-ones(1,6)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])
  x6 = normalize_node(x[6],domain[:,6])

  polynomial_1_data = zeros(2,order+1)
  polynomial_1_data[1,:] = chebyshev_polynomial(order,x1)
  polynomial_1_data[2,:] = chebyshev_polynomial_derivative(order,x1)

  polynomial_2_data = zeros(2,order+1)
  polynomial_2_data[1,:] = chebyshev_polynomial(order,x2)
  polynomial_2_data[2,:] = chebyshev_polynomial_derivative(order,x2)

  polynomial_3_data = zeros(2,order+1)
  polynomial_3_data[1,:] = chebyshev_polynomial(order,x3)
  polynomial_3_data[2,:] = chebyshev_polynomial_derivative(order,x3)

  polynomial_4_data = zeros(2,order+1)
  polynomial_4_data[1,:] = chebyshev_polynomial(order,x4)
  polynomial_4_data[2,:] = chebyshev_polynomial_derivative(order,x4)

  polynomial_5_data = zeros(2,order+1)
  polynomial_5_data[1,:] = chebyshev_polynomial(order,x5)
  polynomial_5_data[2,:] = chebyshev_polynomial_derivative(order,x5)

  polynomial_6_data = zeros(2,order+1)
  polynomial_6_data[1,:] = chebyshev_polynomial(order,x6)
  polynomial_6_data[2,:] = chebyshev_polynomial_derivative(order,x6)

  evaluated_derivative = zeros(1,6)

  for n = 1:order+1
	  for m = 1:order+1
	    for l = 1:order+1
        for k = 1:order+1
          for j = 1:order+1
            for i = 1:order+1

              if (i+j+k+l+m+n <= order+6)

                evaluated_derivative[1] += (2.0/(domain[1,1]-domain[2,1]))*polynomial_1_data[2,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
                evaluated_derivative[2] += (2.0/(domain[1,2]-domain[2,2]))*polynomial_1_data[1,i]*polynomial_2_data[2,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
                evaluated_derivative[3] += (2.0/(domain[1,3]-domain[2,3]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[2,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
                evaluated_derivative[4] += (2.0/(domain[1,4]-domain[2,4]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[2,l]*polynomial_5_data[1,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
                evaluated_derivative[5] += (2.0/(domain[1,5]-domain[2,5]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[2,m]*polynomial_6_data[1,n]*weights[i,j,k,l,m,n]
                evaluated_derivative[6] += (2.0/(domain[1,6]-domain[2,6]))*polynomial_1_data[1,i]*polynomial_2_data[1,j]*polynomial_3_data[1,k]*polynomial_4_data[1,l]*polynomial_5_data[1,m]*polynomial_6_data[2,n]*weights[i,j,k,l,m,n]

              end

            end
          end
				end
      end
    end
  end

  return evaluated_derivative

end

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
"""
Returns a tapered mesh using dDwn and dUp based on the geometry of an input mesh.
Tapering a mesh transforms it into a tapered version (*i.e.* pointy) or enlarged object,
*e.g.* make a cone from a cylinder.
"""
function taper(mesh, dUp, dDwn)
    if dUp != 1.0 && dDwn != 1.0 && !isnan(dUp) && !isnan(dDwn)
        mesh_points = mesh.points
        delta = dDwn - dUp
        Xs = map(x -> x.coords[1], mesh_points)
        xmin = minimum(Xs)
        xmax = maximum(Xs)
        deltaX = xmax - xmin

        scaled_mesh = Array{Meshes.Point3}(undef, length(mesh_points))
        for i = 1:length(mesh_points)
            dX = (mesh_points[i].coords[1] - xmin)
            factor = dDwn - delta * (dX / deltaX)
            scaled_mesh[i] = Meshes.Point3(
                mesh_points[i].coords[1],
                mesh_points[i].coords[2] * factor,
                mesh_points[i].coords[3] * factor
            )
        end
        mesh = Meshes.SimpleMesh(scaled_mesh, Meshes.topology(mesh))
    end

    return mesh
end


text:
@testset "Bondarenko Problems" begin
    (x, z, ) = @variables x z

    #= ::Package:: =#

    #= ::Title:: =#
    #=Vladimir*Bondarenko*Integration*Problems=#

    #= ::Section::Closed:: =#
    #=9*June*2010=#

    @test_int [1/(sqrt(2) + sin(z) + cos(z)), z, 1, -((1 - sqrt(2)*sin(z))/(cos(z) - sin(z)))]

    @test_int [1/(sqrt(1 + x) + sqrt(1 - x))^2, x, 4, -(1/(2*x)) + sqrt(1 - x^2)/(2*x) + asin(x)/2]

    @test_int [1/(1 + cos(x))^2, x, 2, sin(x)/(3*(1 + cos(x))^2) + sin(x)/(3*(1 + cos(x)))]
    @test_int [sin(x)/sqrt(1 + x), x, 5, sqrt(2*pi)*cos(1)*FresnelS(sqrt(2/pi)*sqrt(1 + x)) - sqrt(2*pi)*FresnelC(sqrt(2/pi)*sqrt(1 + x))*sin(1)]
    @test_int [1/(cos(x) + sin(x))^6, x, 3, -((cos(x) - sin(x))/(10*(cos(x) + sin(x))^5)) - (cos(x) - sin(x))/(15*(cos(x) + sin(x))^3) + (2*sin(x))/(15*(cos(x) + sin(x)))]

    @test_int [log(x^4 + 1/x^4), x, 22, -4*x - sqrt(2 + sqrt(2))*atan((sqrt(2 - sqrt(2)) - 2*x)/sqrt(2 + sqrt(2))) - sqrt(2 - sqrt(2))*atan((sqrt(2 + sqrt(2)) - 2*x)/sqrt(2 - sqrt(2))) + sqrt(2 + sqrt(2))*atan((sqrt(2 - sqrt(2)) + 2*x)/sqrt(2 + sqrt(2))) + sqrt(2 - sqrt(2))*atan((sqrt(2 + sqrt(2)) + 2*x)/sqrt(2 - sqrt(2))) - (1/2)*sqrt(2 - sqrt(2))*log(1 - sqrt(2 - sqrt(2))*x + x^2) + (1/2)*sqrt(2 - sqrt(2))*log(1 + sqrt(2 - sqrt(2))*x + x^2) - (1/2)*sqrt(2 + sqrt(2))*log(1 - sqrt(2 + sqrt(2))*x + x^2) + (1/2)*sqrt(2 + sqrt(2))*log(1 + sqrt(2 + sqrt(2))*x + x^2) + x*log(1/x^4 + x^4)]
    @test_int [log(1 + x)/(x*sqrt(1 + sqrt(1 + x))), x, -1, -8*atanh(sqrt(1 + sqrt(1 + x))) - (2*log(1 + x))/sqrt(1 + sqrt(1 + x)) - sqrt(2)*atanh(sqrt(1 + sqrt(1 + x))/sqrt(2))*log(1 + x) + 2*sqrt(2)*atanh(1/sqrt(2))*log(1 - sqrt(1 + sqrt(1 + x))) - 2*sqrt(2)*atanh(1/sqrt(2))*log(1 + sqrt(1 + sqrt(1 + x))) + sqrt(2)*PolyLog(2, -((sqrt(2)*(1 - sqrt(1 + sqrt(1 + x))))/(2 - sqrt(2)))) - sqrt(2)*PolyLog(2, (sqrt(2)*(1 - sqrt(1 + sqrt(1 + x))))/(2 + sqrt(2))) - sqrt(2)*PolyLog(2, -((sqrt(2)*(1 + sqrt(1 + sqrt(1 + x))))/(2 - sqrt(2)))) + sqrt(2)*PolyLog(2, (sqrt(2)*(1 + sqrt(1 + sqrt(1 + x))))/(2 + sqrt(2)))]
    @test_int [log(1 + x)/x*sqrt(1 + sqrt(1 + x)), x, -1, -16*sqrt(1 + sqrt(1 + x)) + 16*atanh(sqrt(1 + sqrt(1 + x))) + 4*sqrt(1 + sqrt(1 + x))*log(1 + x) - 2*sqrt(2)*atanh(sqrt(1 + sqrt(1 + x))/sqrt(2))*log(1 + x) + 4*sqrt(2)*atanh(1/sqrt(2))*log(1 - sqrt(1 + sqrt(1 + x))) - 4*sqrt(2)*atanh(1/sqrt(2))*log(1 + sqrt(1 + sqrt(1 + x))) + 2*sqrt(2)*PolyLog(2, -((sqrt(2)*(1 - sqrt(1 + sqrt(1 + x))))/(2 - sqrt(2)))) - 2*sqrt(2)*PolyLog(2, (sqrt(2)*(1 - sqrt(1 + sqrt(1 + x))))/(2 + sqrt(2))) - 2*sqrt(2)*PolyLog(2, -((sqrt(2)*(1 + sqrt(1 + sqrt(1 + x))))/(2 - sqrt(2)))) + 2*sqrt(2)*PolyLog(2, (sqrt(2)*(1 + sqrt(1 + sqrt(1 + x))))/(2 + sqrt(2)))]

    #= ::Section::Closed:: =#
    #=4*July*2010=#

    @test_int [1/(1 + sqrt(x + sqrt(1 + x^2))), x, 4, -(1/(2*(x + sqrt(1 + x^2)))) + 1/sqrt(x + sqrt(1 + x^2)) + sqrt(x + sqrt(1 + x^2)) + (1/2)*log(x + sqrt(1 + x^2)) - 2*log(1 + sqrt(x + sqrt(1 + x^2)))]
    @test_int [sqrt(1 + x)/(x + sqrt(1 + sqrt(1 + x))), x, 6, 2*sqrt(1 + x) + (8*atanh((1 + 2*sqrt(1 + sqrt(1 + x)))/sqrt(5)))/sqrt(5)]
    @test_int [1/(x - sqrt(1 + sqrt(1 + x))), x, 5, (2/5)*(5 + sqrt(5))*log(1 - sqrt(5) - 2*sqrt(1 + sqrt(1 + x))) + (2/5)*(5 - sqrt(5))*log(1 + sqrt(5) - 2*sqrt(1 + sqrt(1 + x)))]
    @test_int [x/(x + sqrt(1 - sqrt(1 + x))), x, 6, 2*sqrt(1 + x) - 4*sqrt(1 - sqrt(1 + x)) + (1 - sqrt(1 + x))^2 + (8*atanh((1 + 2*sqrt(1 - sqrt(1 + x)))/sqrt(5)))/sqrt(5)]
    @test_int [sqrt(sqrt(1 + x) + x)/((1 + x^2)*sqrt(1 + x)), x, 20, -((im*atan((2 + sqrt(1 - im) - (1 - 2*sqrt(1 - im))*sqrt(1 + x))/(2*sqrt(im + sqrt(1 - im))*sqrt(x + sqrt(1 + x)))))/(2*sqrt((1 - im)/(im + sqrt(1 - im))))) + (im*atan((2 + sqrt(1 + im) - (1 - 2*sqrt(1 + im))*sqrt(1 + x))/(2*sqrt(-im + sqrt(1 + im))*sqrt(x + sqrt(1 + x)))))/(2*sqrt(-((1 + im)/(im - sqrt(1 + im))))) + (im*atanh((2 - sqrt(1 - im) - (1 + 2*sqrt(1 - im))*sqrt(1 + x))/(2*sqrt(-im + sqrt(1 - im))*sqrt(x + sqrt(1 + x)))))/(2*sqrt(-((1 - im)/(im - sqrt(1 - im))))) - (im*atanh((2 - sqrt(1 + im) - (1 + 2*sqrt(1 + im))*sqrt(1 + x))/(2*sqrt(im + sqrt(1 + im))*sqrt(x + sqrt(1 + x)))))/(2*sqrt((1 + im)/(im + sqrt(1 + im))))]
    @test_int [sqrt(x + sqrt(1 + x))/(1 + x^2), x, 22, (1/2)*im*sqrt(im + sqrt(1 - im))*atan((2 + sqrt(1 - im) - (1 - 2*sqrt(1 - im))*sqrt(1 + x))/(2*sqrt(im + sqrt(1 - im))*sqrt(x + sqrt(1 + x)))) - (1/2)*im*sqrt(-im + sqrt(1 + im))*atan((2 + sqrt(1 + im) - (1 - 2*sqrt(1 + im))*sqrt(1 + x))/(2*sqrt(-im + sqrt(1 + im))*sqrt(x + sqrt(1 + x)))) + (1/2)*im*sqrt(-im + sqrt(1 - im))*atanh((2 - sqrt(1 - im) - (1 + 2*sqrt(1 - im))*sqrt(1 + x))/(2*sqrt(-im + sqrt(1 - im))*sqrt(x + sqrt(1 + x)))) - (1/2)*im*sqrt(im + sqrt(1 + im))*atanh((2 - sqrt(1 + im) - (1 + 2*sqrt(1 + im))*sqrt(1 + x))/(2*sqrt(im + sqrt(1 + im))*sqrt(x + sqrt(1 + x))))]
    @test_int [sqrt(1 + sqrt(x) + sqrt(1 + 2*sqrt(x) + 2*x)), x, 2, (2*sqrt(1 + sqrt(x) + sqrt(1 + 2*sqrt(x) + 2*x))*(2 + sqrt(x) + 6*x^(3/2) - (2 - sqrt(x))*sqrt(1 + 2*sqrt(x) + 2*x)))/(15*sqrt(x))]
    @test_int [sqrt(sqrt(2) + sqrt(x) + sqrt(2 + sqrt(8)*sqrt(x) + 2*x)), x, 3, (1/(15*sqrt(x)))*(2*sqrt(2)*sqrt(sqrt(2) + sqrt(x) + sqrt(2)*sqrt(1 + sqrt(2)*sqrt(x) + x))*(4 + sqrt(2)*sqrt(x) + 3*sqrt(2)*x^(3/2) - sqrt(2)*(2*sqrt(2) - sqrt(x))*sqrt(1 + sqrt(2)*sqrt(x) + x)))]
    @test_int [sqrt(x + sqrt(1 + x))/x^2, x, 7, -(sqrt(x + sqrt(1 + x))/x) - (1/4)*atan((3 + sqrt(1 + x))/(2*sqrt(x + sqrt(1 + x)))) + (3/4)*atanh((1 - 3*sqrt(1 + x))/(2*sqrt(x + sqrt(1 + x))))]
    @test_int [sqrt(1/x + sqrt(1 + 1/x)), x, 7, sqrt(sqrt(1 + 1/x) + 1/x)*x + (1/4)*atan((3 + sqrt(1 + 1/x))/(2*sqrt(sqrt(1 + 1/x) + 1/x))) - (3/4)*atanh((1 - 3*sqrt(1 + 1/x))/(2*sqrt(sqrt(1 + 1/x) + 1/x)))]

    @test_int [sqrt(1 + exp(-x))/(exp(x) - exp(-x)), x, 6, (-sqrt(2))*atanh(sqrt(1 + ℯ^(-x))/sqrt(2))]
    @test_int [sqrt(1 + exp(-x))/sinh(x), x, 7, -2*sqrt(2)*atanh(sqrt(1 + ℯ^(-x))/sqrt(2))]

    @test_int [1/(cos(x) + cos(3*x))^5, x, -45, (-(523/256))*atanh(sin(x)) + (1483*atanh(sqrt(2)*sin(x)))/(512*sqrt(2)) + sin(x)/(32*(1 - 2*sin(x)^2)^4) - (17*sin(x))/(192*(1 - 2*sin(x)^2)^3) + (203*sin(x))/(768*(1 - 2*sin(x)^2)^2) - (437*sin(x))/(512*(1 - 2*sin(x)^2)) - (43/256)*sec(x)*tan(x) - (1/128)*sec(x)^3*tan(x)]
    @test_int [1/(cos(x) + sin(x) + 1)^2, x, 3, -log(1 + tan(x/2)) - (cos(x) - sin(x))/(1 + cos(x) + sin(x))]

    @test_int [sqrt(1 + tanh(4*x)), x, 2, atanh(sqrt(1 + tanh(4*x))/sqrt(2))/(2*sqrt(2))]
    @test_int [tanh(x)/sqrt(exp(2*x) + exp(x)), x, -11, (2*sqrt(ℯ^x + ℯ^(2*x)))/ℯ^x - atan((im - (1 - 2*im)*ℯ^x)/(2*sqrt(1 + im)*sqrt(ℯ^x + ℯ^(2*x))))/sqrt(1 + im) + atan((im + (1 + 2*im)*ℯ^x)/(2*sqrt(1 - im)*sqrt(ℯ^x + ℯ^(2*x))))/sqrt(1 - im)]
    @test_int [sqrt(sinh(2*x)/cosh(x)), x, 5, (2*im*sqrt(2)*Elliptic.E(pi/4 - (im*x)/2, 2)*sqrt(sinh(x)))/sqrt(im*sinh(x)), (2*im*Elliptic.E(pi/4 - (im*x)/2, 2)*sqrt(sech(x)*sinh(2*x)))/sqrt(im*sinh(x))]

    @test_int [log(x^2 + sqrt(1 - x^2)), x, -31, -2*x - asin(x) + sqrt((1/2)*(1 + sqrt(5)))*atan(sqrt(2/(1 + sqrt(5)))*x) + sqrt((1/2)*(1 + sqrt(5)))*atan((sqrt((1/2)*(1 + sqrt(5)))*x)/sqrt(1 - x^2)) + sqrt((1/2)*(-1 + sqrt(5)))*atanh(sqrt(2/(-1 + sqrt(5)))*x) - sqrt((1/2)*(-1 + sqrt(5)))*atanh((sqrt((1/2)*(-1 + sqrt(5)))*x)/sqrt(1 - x^2)) + x*log(x^2 + sqrt(1 - x^2))]
    @test_int [log(1 + exp(x))/(1 + exp(2*x)), x, 12, (-(1/2))*log((1/2 - im/2)*(im - ℯ^x))*log(1 + ℯ^x) - (1/2)*log((-(1/2) - im/2)*(im + ℯ^x))*log(1 + ℯ^x) - PolyLog(2, -ℯ^x) - (1/2)*PolyLog(2, (1/2 - im/2)*(1 + ℯ^x)) - (1/2)*PolyLog(2, (1/2 + im/2)*(1 + ℯ^x))]
    @test_int [log(1 + cosh(x)^2)^2*cosh(x), x, 13, -8*sqrt(2)*atan(sinh(x)/sqrt(2)) + 4*im*sqrt(2)*atan(sinh(x)/sqrt(2))^2 + 8*sqrt(2)*atan(sinh(x)/sqrt(2))*log((2*sqrt(2))/(sqrt(2) + im*sinh(x))) + 4*sqrt(2)*atan(sinh(x)/sqrt(2))*log(2 + sinh(x)^2) + 4*im*sqrt(2)*PolyLog(2, 1 - (2*sqrt(2))/(sqrt(2) + im*sinh(x))) + 8*sinh(x) - 4*log(2 + sinh(x)^2)*sinh(x) + log(2 + sinh(x)^2)^2*sinh(x)]
    @test_int [log(sinh(x) + cosh(x)^2)^2*cosh(x), x, 28, -4*sqrt(3)*atan((1 + 2*sinh(x))/sqrt(3)) - (1/2)*(1 - im*sqrt(3))*log(1 - im*sqrt(3) + 2*sinh(x))^2 - (1 + im*sqrt(3))*log((im*(1 - im*sqrt(3) + 2*sinh(x)))/(2*sqrt(3)))*log(1 + im*sqrt(3) + 2*sinh(x)) - (1/2)*(1 + im*sqrt(3))*log(1 + im*sqrt(3) + 2*sinh(x))^2 - (1 - im*sqrt(3))*log(1 - im*sqrt(3) + 2*sinh(x))*log(-((im*(1 + im*sqrt(3) + 2*sinh(x)))/(2*sqrt(3)))) - 2*log(1 + sinh(x) + sinh(x)^2) + (1 - im*sqrt(3))*log(1 - im*sqrt(3) + 2*sinh(x))*log(1 + sinh(x) + sinh(x)^2) + (1 + im*sqrt(3))*log(1 + im*sqrt(3) + 2*sinh(x))*log(1 + sinh(x) + sinh(x)^2) - (1 + im*sqrt(3))*PolyLog(2, -((im - sqrt(3) + 2*im*sinh(x))/(2*sqrt(3)))) - (1 - im*sqrt(3))*PolyLog(2, (im + sqrt(3) + 2*im*sinh(x))/(2*sqrt(3))) + 8*sinh(x) - 4*log(1 + sinh(x) + sinh(x)^2)*sinh(x) + log(1 + sinh(x) + sinh(x)^2)^2*sinh(x)]
    @test_int [log(x + sqrt(1 + x))/(1 + x^2), x, 44, (1/2)*im*log(sqrt(1 - im) - sqrt(1 + x))*log(x + sqrt(1 + x)) - (1/2)*im*log(sqrt(1 + im) - sqrt(1 + x))*log(x + sqrt(1 + x)) + (1/2)*im*log(sqrt(1 - im) + sqrt(1 + x))*log(x + sqrt(1 + x)) - (1/2)*im*log(sqrt(1 + im) + sqrt(1 + x))*log(x + sqrt(1 + x)) - (1/2)*im*log(sqrt(1 - im) + sqrt(1 + x))*log((1 - sqrt(5) + 2*sqrt(1 + x))/(1 - 2*sqrt(1 - im) - sqrt(5))) - (1/2)*im*log(sqrt(1 - im) - sqrt(1 + x))*log((1 - sqrt(5) + 2*sqrt(1 + x))/(1 + 2*sqrt(1 - im) - sqrt(5))) + (1/2)*im*log(sqrt(1 + im) + sqrt(1 + x))*log((1 - sqrt(5) + 2*sqrt(1 + x))/(1 - 2*sqrt(1 + im) - sqrt(5))) + (1/2)*im*log(sqrt(1 + im) - sqrt(1 + x))*log((1 - sqrt(5) + 2*sqrt(1 + x))/(1 + 2*sqrt(1 + im) - sqrt(5))) - (1/2)*im*log(sqrt(1 - im) + sqrt(1 + x))*log((1 + sqrt(5) + 2*sqrt(1 + x))/(1 - 2*sqrt(1 - im) + sqrt(5))) - (1/2)*im*log(sqrt(1 - im) - sqrt(1 + x))*log((1 + sqrt(5) + 2*sqrt(1 + x))/(1 + 2*sqrt(1 - im) + sqrt(5))) + (1/2)*im*log(sqrt(1 + im) + sqrt(1 + x))*log((1 + sqrt(5) + 2*sqrt(1 + x))/(1 - 2*sqrt(1 + im) + sqrt(5))) + (1/2)*im*log(sqrt(1 + im) - sqrt(1 + x))*log((1 + sqrt(5) + 2*sqrt(1 + x))/(1 + 2*sqrt(1 + im) + sqrt(5))) - (1/2)*im*PolyLog(2, (2*(sqrt(1 - im) - sqrt(1 + x)))/(1 + 2*sqrt(1 - im) - sqrt(5))) - (1/2)*im*PolyLog(2, (2*(sqrt(1 - im) - sqrt(1 + x)))/(1 + 2*sqrt(1 - im) + sqrt(5))) + (1/2)*im*PolyLog(2, (2*(sqrt(1 + im) - sqrt(1 + x)))/(1 + 2*sqrt(1 + im) - sqrt(5))) + (1/2)*im*PolyLog(2, (2*(sqrt(1 + im) - sqrt(1 + x)))/(1 + 2*sqrt(1 + im) + sqrt(5))) - (1/2)*im*PolyLog(2, -((2*(sqrt(1 - im) + sqrt(1 + x)))/(1 - 2*sqrt(1 - im) - sqrt(5)))) - (1/2)*im*PolyLog(2, -((2*(sqrt(1 - im) + sqrt(1 + x)))/(1 - 2*sqrt(1 - im) + sqrt(5)))) + (1/2)*im*PolyLog(2, -((2*(sqrt(1 + im) + sqrt(1 + x)))/(1 - 2*sqrt(1 + im) - sqrt(5)))) + (1/2)*im*PolyLog(2, -((2*(sqrt(1 + im) + sqrt(1 + x)))/(1 - 2*sqrt(1 + im) + sqrt(5))))]
    @test_int [log(x + sqrt(1 + x))^2/(1 + x)^2, x, 35, log(1 + x) + (2*log(x + sqrt(1 + x)))/sqrt(1 + x) - 6*log(sqrt(1 + x))*log(x + sqrt(1 + x)) - log(x + sqrt(1 + x))^2/(1 + x) - (1 + sqrt(5))*log(1 - sqrt(5) + 2*sqrt(1 + x)) + 6*log((1/2)*(-1 + sqrt(5)))*log(1 - sqrt(5) + 2*sqrt(1 + x)) + (3 + sqrt(5))*log(x + sqrt(1 + x))*log(1 - sqrt(5) + 2*sqrt(1 + x)) - (1/2)*(3 + sqrt(5))*log(1 - sqrt(5) + 2*sqrt(1 + x))^2 - (1 - sqrt(5))*log(1 + sqrt(5) + 2*sqrt(1 + x)) + (3 - sqrt(5))*log(x + sqrt(1 + x))*log(1 + sqrt(5) + 2*sqrt(1 + x)) - (3 - sqrt(5))*log(-((1 - sqrt(5) + 2*sqrt(1 + x))/(2*sqrt(5))))*log(1 + sqrt(5) + 2*sqrt(1 + x)) - (1/2)*(3 - sqrt(5))*log(1 + sqrt(5) + 2*sqrt(1 + x))^2 - (3 + sqrt(5))*log(1 - sqrt(5) + 2*sqrt(1 + x))*log((1 + sqrt(5) + 2*sqrt(1 + x))/(2*sqrt(5))) + 6*log(sqrt(1 + x))*log(1 + (2*sqrt(1 + x))/(1 + sqrt(5))) + 6*PolyLog(2, -((2*sqrt(1 + x))/(1 + sqrt(5)))) - (3 + sqrt(5))*PolyLog(2, -((1 - sqrt(5) + 2*sqrt(1 + x))/(2*sqrt(5)))) - (3 - sqrt(5))*PolyLog(2, (1 + sqrt(5) + 2*sqrt(1 + x))/(2*sqrt(5))) - 6*PolyLog(2, 1 + (2*sqrt(1 + x))/(1 - sqrt(5)))]
    @test_int [log(x + sqrt(1 + x))/x, x, 21, log(-1 + sqrt(1 + x))*log(x + sqrt(1 + x)) + log(1 + sqrt(1 + x))*log(x + sqrt(1 + x)) - log(-1 + sqrt(1 + x))*log((1 - sqrt(5) + 2*sqrt(1 + x))/(3 - sqrt(5))) - log(1 + sqrt(1 + x))*log(-((1 - sqrt(5) + 2*sqrt(1 + x))/(1 + sqrt(5)))) - log(1 + sqrt(1 + x))*log(-((1 + sqrt(5) + 2*sqrt(1 + x))/(1 - sqrt(5)))) - log(-1 + sqrt(1 + x))*log((1 + sqrt(5) + 2*sqrt(1 + x))/(3 + sqrt(5))) - PolyLog(2, (2*(1 - sqrt(1 + x)))/(3 - sqrt(5))) - PolyLog(2, (2*(1 - sqrt(1 + x)))/(3 + sqrt(5))) - PolyLog(2, (2*(1 + sqrt(1 + x)))/(1 - sqrt(5))) - PolyLog(2, (2*(1 + sqrt(1 + x)))/(1 + sqrt(5)))]

    @test_int [atan(2*tan(x)), x, 7, x*atan(2*tan(x)) + (1/2)*im*x*log(1 - 3*ℯ^(2*im*x)) - (1/2)*im*x*log(1 - (1/3)*ℯ^(2*im*x)) - (1/4)*PolyLog(2, (1/3)*ℯ^(2*im*x)) + (1/4)*PolyLog(2, 3*ℯ^(2*im*x))]
    @test_int [atan(x)*log(x)/x, x, 5, (1/2)*im*log(x)*PolyLog(2, (-im)*x) - (1/2)*im*log(x)*PolyLog(2, im*x) - (1/2)*im*PolyLog(3, (-im)*x) + (1/2)*im*PolyLog(3, im*x)]

    #= Note: Mathematica*is*unable*to*differentiate*result*back*to*integrand! =#
    @test_int [atan(x)^2*sqrt(1 + x^2), x, 10, asinh(x) - sqrt(1 + x^2)*atan(x) + (1/2)*x*sqrt(1 + x^2)*atan(x)^2 - im*atan(ℯ^(im*atan(x)))*atan(x)^2 + im*atan(x)*PolyLog(2, (-im)*ℯ^(im*atan(x))) - im*atan(x)*PolyLog(2, im*ℯ^(im*atan(x))) - PolyLog(3, (-im)*ℯ^(im*atan(x))) + PolyLog(3, im*ℯ^(im*atan(x)))]
    end

max_stars_repo_path: test/0 Independent test suites/Bondarenko Problems.jl


text:
"""
    flux_variability_analysis(
        model::MetabolicModel,
        reactions::Vector{Int},
        optimizer;
        modifications = [],
        workers = [myid()],
        optimal_objective_value = nothing,
        bounds = z -> (z, Inf),
        ret = objective_value,
    )::Matrix{Float64}

Flux variability analysis solves a pair of optimization problems in `model` for
each flux listed in `reactions`:
```
 min,max xᵢ
s.t. S x = b
    xₗ ≤ x ≤ xᵤ
     cᵀx ≥ bounds(Z₀)[1]
     cᵀx ≤ bounds(Z₀)[2]
```
where Z₀:= cᵀx₀ is the objective value of an optimal solution of the associated
FBA problem (see [`flux_balance_analysis`](@ref) for a related analysis, also
for explanation of the `modifications` argument).

The `bounds` is a user-supplied function that specifies the objective bounds
for the variability optimizations, by default it restricts the flux objective
value to the precise optimum reached in FBA. It can return `-Inf` and `Inf` in
first and second pair to remove the limit. Use [`gamma_bounds`](@ref) and
[`objective_bounds`](@ref) for simple bounds.

`optimizer` must be set to a `JuMP`-compatible optimizer. The computation of
the individual optimization problems is transparently distributed to `workers`
(see `Distributed.workers()`).  The value of Z₀ can be optionally supplied in
argument `optimal_objective_value`, which prevents this function from calling
the non-parallelizable FBA. Separating the single-threaded FBA and
multithreaded variability computation can be used to improve resource
allocation efficiency in many common use-cases.

`ret` is a function used to extract results from optimized JuMP models of the
individual reactions. By default, it calls and returns the value of
`JuMP.objective_value`. More information can be extracted e.g. by setting it to
a function that returns a more elaborate data structure; such as `m ->
(JuMP.objective_value(m), JuMP.value.(m[:x]))`.

Returns a matrix of extracted `ret` values for minima and maxima, of total size
(`length(reactions)`,2). The optimizer result status is checked with
[`is_solved`](@ref); `nothing` is returned if the optimization failed for any
reason.

# Example
```
model = load_model("e_coli_core.json")
flux_variability_analysis(model, [1, 2, 3, 42], GLPK.optimizer)
```
"""
function flux_variability_analysis(
    model::MetabolicModel,
    reactions::Vector{Int},
    optimizer;
    modifications = [],
    workers = [myid()],
    optimal_objective_value = nothing,
    bounds = z -> (z, Inf),
    ret = objective_value,
)
    if any(reactions .< 1) || any(reactions .> n_reactions(model))
        throw(DomainError(reactions, "Index exceeds number of reactions."))
    end

    Z = bounds(
        isnothing(optimal_objective_value) ?
        objective_value(
            flux_balance_analysis(model, optimizer; modifications = modifications),
        ) : optimal_objective_value,
    )

    return screen_optmodel_modifications(
        model,
        optimizer;
        common_modifications = vcat(
            modifications,
            [
                (model, opt_model) -> begin
                    Z[1] > -Inf && @constraint(
                        opt_model,
                        objective(model)' * opt_model[:x] >= Z[1]
                    )
                    Z[2] < Inf && @constraint(
                        opt_model,
                        objective(model)' * opt_model[:x] <= Z[2]
                    )
                end,
            ],
        ),
        args = Tuple.([-reactions reactions]),
        analysis = (_, opt_model, ridx) -> _max_variability_flux(opt_model, ridx, ret),
        workers = workers,
    )
end

"""
    flux_variability_analysis(
        model::MetabolicModel,
        optimizer;
        kwargs...
    )

A simpler version of [`flux_variability_analysis`](@ref) that maximizes and
minimizes all reactions in the model. Arguments are forwarded.
"""
function flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...)
    n = n_reactions(model)
    return flux_variability_analysis(model, collect(1:n), optimizer; kwargs...)
end

"""
    flux_variability_analysis_dict(
        model::MetabolicModel,
        optimizer;
        kwargs...
    )

A variant of [`flux_variability_analysis`](@ref) that returns the individual
maximized and minimized fluxes of all reactions as two dictionaries (of
dictionaries). All keyword arguments except `ret` are passed through.

# Example
```
mins, maxs = flux_variability_analysis_dict(
    model,
    Tulip.Optimizer;
    bounds = objective_bounds(0.99),
    modifications = [
        change_optimizer_attribute("IPM_IterationsLimit", 500),
        change_constraint("EX_glc__D_e"; lb = -10, ub = -10),
        change_constraint("EX_o2_e"; lb = 0, ub = 0),
    ],
)
```
"""
function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
    fluxes = flux_variability_analysis(model, optimizer; kwargs..., ret = flux_vector)
    rxns = reactions(model)
    dicts = zip.(Ref(rxns), fluxes)

    return (Dict(rxns .=> Dict.(dicts[:, 1])), Dict(rxns .=> Dict.(dicts[:, 2])))
end

"""
    _max_variability_flux(opt_model, rid, ret)

Internal helper for maximizing reactions in optimization model.
"""
function _max_variability_flux(opt_model, rid, ret)
    sense = rid > 0 ? MAX_SENSE : MIN_SENSE
    var = all_variables(opt_model)[abs(rid)]

    @objective(opt_model, sense, var)
    optimize!(opt_model)

    is_solved(opt_model) ? ret(opt_model) : nothing
end


text:
vec2mat(v::AbstractVector) = reshape(v, (length(v),1))
vec2mat(T::Type, v::AbstractVector) = vec2mat(T.(v))