text:
# scalar version

"""
    ScalarSlackBridge{T, F, S}

The `ScalarSlackBridge` converts a constraint `G`-in-`S` where `G` is a function different
from `SingleVariable` into the constraints `F`-in-`EqualTo{T}` and `SingleVariable`-in-`S`.
`F` is the result of subtracting a `SingleVariable` from `G`.
Tipically `G` is the same as `F`, but that is not mandatory.
"""
struct ScalarSlackBridge{T, F, S} <: AbstractBridge
    slack::MOI.VariableIndex
    slack_in_set::CI{MOI.SingleVariable, S}
    equality::CI{F, MOI.EqualTo{T}}
end
function ScalarSlackBridge{T, F, S}(model, f::MOI.AbstractScalarFunction, s::S) where {T, F, S}
    slack = MOI.add_variable(model)
    new_f = MOIU.operate(-, T, f, MOI.SingleVariable(slack))
    slack_in_set = MOI.add_constraint(model, MOI.SingleVariable(slack), s)
    equality = MOI.add_constraint(model, new_f, MOI.EqualTo(0.0))
    return ScalarSlackBridge{T, F, S}(slack, slack_in_set, equality)
end

# start allowing everything (scalar)
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.AbstractScalarFunction},
                        ::Type{<:MOI.AbstractScalarSet}) where {T} = true
# then restrict (careful with ambiguity)
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.SingleVariable},
                        ::Type{<:MOI.EqualTo}) where {T} = false
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.SingleVariable},
                        ::Type{<:MOI.AbstractScalarSet}) where {T} = false
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.AbstractScalarFunction},
                        ::Type{<:MOI.EqualTo}) where {T} = false
function added_constraint_types(::Type{ScalarSlackBridge{T, F, S}}) where {T, F, S}
    return [(F, MOI.EqualTo{T}), (MOI.SingleVariable, S)]
end
function concrete_bridge_type(::Type{<:ScalarSlackBridge{T}},
                              F::Type{<:MOI.AbstractScalarFunction},
                              S::Type{<:MOI.AbstractScalarSet}) where T
    F2 = MOIU.promote_operation(-, T, F, MOI.SingleVariable)
    return ScalarSlackBridge{T, F2, S}
end

# Attributes, Bridge acting as an model
MOI.get(b::ScalarSlackBridge{T, F, S}, ::MOI.NumberOfVariables) where {T, F, S} = 1
MOI.get(b::ScalarSlackBridge{T, F, S}, ::MOI.NumberOfConstraints{F, MOI.EqualTo{T}}) where {T, F, S} = 1
MOI.get(b::ScalarSlackBridge{T, F, S}, ::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where {T, F, S} = 1
MOI.get(b::ScalarSlackBridge{T, F, S}, ::MOI.ListOfConstraintIndices{F, MOI.EqualTo{T}}) where {T, F, S} = [b.equality]
MOI.get(b::ScalarSlackBridge{T, F, S}, ::MOI.ListOfConstraintIndices{MOI.SingleVariable, S}) where {T, F, S} = [b.slack_in_set]

# Indices
function MOI.delete(model::MOI.ModelLike, c::ScalarSlackBridge)
    MOI.delete(model, c.equality)
    MOI.delete(model, c.slack_in_set)
    MOI.delete(model, c.slack)
    return
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal, c::ScalarSlackBridge)
    # due to equality, slack should have the same value as original affine function
    return MOI.get(model, attr, c.slack_in_set)
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintDual, c::ScalarSlackBridge)
    # The dual constraint on slack (since it is free) is
    # -dual_slack_in_set + dual_equality = 0 so the two duals are
    # equal and we can return either one of them.
    return MOI.get(model, a, c.slack_in_set)
end
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintBasisStatus,  c::ScalarSlackBridge)
    MOI.get(model, MOI.ConstraintBasisStatus(), c.slack_in_set)
end

# Constraints
function MOI.modify(model::MOI.ModelLike, c::ScalarSlackBridge, change::MOI.AbstractFunctionModification)
    MOI.modify(model, c.equality, change)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 c::ScalarSlackBridge{T, F, S}, func::F) where {T, F, S}
    new_func = MOIU.operate(-, T, func, MOI.SingleVariable(c.slack))
    MOI.set(model, MOI.ConstraintFunction(), c.equality, new_func)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet, c::ScalarSlackBridge{T, F, S}, change::S) where {T, F, S}
    MOI.set(model, MOI.ConstraintSet(), c.slack_in_set, change)
end

# vector version

"""
    VectorSlackBridge{T, F, S}

The `VectorSlackBridge` converts a constraint `G`-in-`S` where `G` is a function different
from `VectorOfVariables` into the constraints `F`in-`Zeros` and `VectorOfVariables`-in-`S`.
`F` is the result of subtracting a `VectorOfVariables` from `G`.
Tipically `G` is the same as `F`, but that is not mandatory.
"""
struct VectorSlackBridge{T, F, S} <: AbstractBridge
    slacks::Vector{MOI.VariableIndex}
    slacks_in_set::CI{MOI.VectorOfVariables, S}
    equality::CI{F, MOI.Zeros}
end
function VectorSlackBridge{T, F, S}(model, f::MOI.AbstractVectorFunction, s::S) where {T, F, S}
    d = MOI.dimension(s)
    slacks = MOI.add_variables(model, d)
    new_f = MOIU.operate(-, T, f, MOI.VectorAffineFunction{T}(MOI.VectorOfVariables(slacks)))
    slacks_in_set = MOI.add_constraint(model, MOI.VectorOfVariables(slacks), s)
    equality = MOI.add_constraint(model, new_f, MOI.Zeros(d))
    return VectorSlackBridge{T, F, S}(slacks, slacks_in_set, equality)
end

MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.AbstractVectorFunction},
                        ::Type{<:MOI.AbstractVectorSet}) where {T} = true
MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.VectorOfVariables},
                        ::Type{<:MOI.Zeros}) where {T} = false
MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.AbstractVectorFunction},
                        ::Type{<:MOI.Zeros}) where {T} = false
MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.VectorOfVariables},
                        ::Type{<:MOI.AbstractVectorSet}) where {T} = false
function added_constraint_types(::Type{VectorSlackBridge{T, F, S}}) where {T, F<:MOI.AbstractVectorFunction, S}
    return [(F, MOI.Zeros), (MOI.VectorOfVariables, S)]
end
function concrete_bridge_type(::Type{<:VectorSlackBridge{T}},
                              F::Type{<:MOI.AbstractVectorFunction},
                              S::Type{<:MOI.AbstractVectorSet}) where T
    F2 = MOIU.promote_operation(-, T, F, MOI.VectorOfVariables)
    return VectorSlackBridge{T, F2, S}
end

# Attributes, Bridge acting as an model
MOI.get(b::VectorSlackBridge{T, F, S}, ::MOI.NumberOfVariables) where {T, F, S} = length(b.slacks)
MOI.get(b::VectorSlackBridge{T, F, S}, ::MOI.NumberOfConstraints{F, MOI.Zeros}) where {T, F, S} = 1
MOI.get(b::VectorSlackBridge{T, F, S}, ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T, F, S} = 1
MOI.get(b::VectorSlackBridge{T, F, S}, ::MOI.ListOfConstraintIndices{F, MOI.Zeros}) where {T, F, S} = [b.equality]
MOI.get(b::VectorSlackBridge{T, F, S}, ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T, F, S} = [b.slacks_in_set]

# Indices
function MOI.delete(model::MOI.ModelLike, c::VectorSlackBridge)
    MOI.delete(model, c.equality)
    MOI.delete(model, c.slacks_in_set)
    MOI.delete(model, c.slacks)
    return
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal, c::VectorSlackBridge)
    # due to equality, slacks should have the same value as original affine function
    return MOI.get(model, attr, c.slacks_in_set)
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintDual, c::VectorSlackBridge)
    # The dual constraint on slack (since it is free) is
    # -dual_slack_in_set + dual_equality = 0 so the two duals are
    # equal and we can return either one of them.
    return MOI.get(model, a, c.slacks_in_set)
end

# Constraints
function MOI.modify(model::MOI.ModelLike, c::VectorSlackBridge, change::MOI.AbstractFunctionModification)
    MOI.modify(model, c.equality, change)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 c::VectorSlackBridge{T, F, S}, func::F) where {T, F, S}
    new_func = MOIU.operate(-, T, func, MOI.VectorAffineFunction{T}(MOI.VectorOfVariables(c.slacks)))
    MOI.set(model, MOI.ConstraintFunction(), c.equality, new_func)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet, c::VectorSlackBridge{T,F,S}, change::S)  where {T, F, S}
    MOI.set(model, MOI.ConstraintSet(), c.slacks_in_set, change)
end


text:
#=
"""
    Miranda(h::Function,X::Vector{Interval{T}},P::Vector{Interval{T}})

Runs Miranda's tests for exclusion of implicit function. Implicit function
defined by `h(X,P) = 0` in (X,P) implies 0 in `{h(x,p} for x in X, p in P}. Inputs:
* `h::Function - Equations defining potential implicit function
* `X::Vector{Interval{T}}` - Bounds for dependent variables
* `P::Vector{Interval{T}}` - Bounds for independent variables
Returns -1 if no exclusion found. Returns index at which exclusion fails otherwise.
"""
function Miranda(h::Function,X::Vector{Interval{T}},P::Vector{Interval{T}}) where {T}
  H:Vector{Interval{T}} = h(X,P)
  for i=1:length(X)
    if (H[i].lo*H[i].hi>0.0)
      return i
    end
  end
  i = -1
  return i
end
=#
#=
"""
    MirandaExc

Runs Miranda's tests for exclusion of implicit function. Implicit function
defined by h(X,P) = 0 in (X,P) implies 0 in {h(x,p} for x in X, p in P}.

Inputs:
* h:         Function - Equations defining potential implicit function
* X:         Vector{Interval{Float64}} - Bounds for dependent variables
* P:         Vector{Interval{Float64}} - Bounds for independent variables
* Eflag_in:  Bool - Input exclusion flag.

Returns true if exclusion found. Returns input otherwise.
"""
function MirandaExc(h,X,P,Eflag_in)
   Eflag = Eflag_in
   i = Miranda(h,X,P)
   if (i != -1)
     Eflag = true
   end
   return Eflag
end
=#
#=
"""
--------------------------------------------------------------------------------
Function: partialIncTop
--------------------------------------------------------------------------------
Description:
Sets intervalbox to an the degenerate interval containing only the upper bound
for dimension i and then apply Miranda's test. Advance to next dimension and
repeat
--------------------------------------------------------------------------------
Inputs:
* h:         function - Equations defining potential implicit function
* X:         Vector{Interval{Float64}} - Bounds for dependent variables
* P:         Vector{Interval{Float64}} - Bounds for independent variables
* PIflag:    Bool - Input exclusion flag.
* incHigh:   Array{Bool} - Array that is false if potential partial inclusion
--------------------------------------------------------------------------------
Returns:
Returns tuple (k,PIflag):
* k:        Int64 - Dimension at which partial inclusion potential exists
* PIflag:   Bool - Partial inclusion flag (false if no partial inclusion)
--------------------------------------------------------------------------------
"""
function partialIncTop(h,X,P,PIflag,incHigh)
  xU = [X[i].hi for i=1:length(X)]
  XU = copy(X)
  Xprev = copy(X)
  k = -1
  ku = -1
  for i=1:length(X)
    XU[i] = Interval(xU[i])
    ku = Miranda(h,XU,P)
    if (ku == -1 && ~incHigh[i])
      k = i
      break
    else
      PIflag = false
    end
    XU[i] = Xprev[i]
  end
  return k,PIflag
end
=#
#=
"""
--------------------------------------------------------------------------------
Function: partialIncBot
--------------------------------------------------------------------------------
Description:
Sets intervalbox to an the degenerate interval containing only the lower bound
for dimension i and then apply Miranda's test. Advance to next dimension and
repeat
--------------------------------------------------------------------------------
Inputs:
* h:         function - Equations defining potential implicit function
* X:         Vector{Interval{Float64}} - Bounds for dependent variables
* P:         Vector{Interval{Float64}} - Bounds for independent variables
* PIflag:    Bool - Input exclusion flag.
* incLow:    Vector{Bool} - Array that is false if potential partial inclusion
--------------------------------------------------------------------------------
Returns:
Returns tuple (k,PIflag):
* k:        Dimension at which partial inclusion potential exists
* PIflag:   Partial inclusion flag (false if no partial inclusion)
--------------------------------------------------------------------------------
"""
function partialIncBot(h,X,P,PIflag,incLow)
  xL = [X[i].lo for i=1:length(X)]
  XL = copy(X)
  Xprev = copy(X)
  k = -1
  kl = -1
  for i=1:length(X)
    XL[i] = Interval(xL[i])
    kl = Miranda(h,XL,P)
    if (kl == -1 && ~incLow[i])
        k = i
      break
    else
      PIflag = false
    end
    XL[i] = Xprev[i]
  end
  return k,PIflag
end
=#
#=
"""
--------------------------------------------------------------------------------
Function: spectralR
--------------------------------------------------------------------------------
Description:
Computes spectral radius for novel inclusion test.
--------------------------------------------------------------------------------
Inputs:
* hj:        function - Jacobian wrt x of potential implicit function
* X:         Vector{Interval{Float64}} - Bounds for dependent variables
* P:         Vector{Interval{Float64}} - Bounds for independent variables
* rho:       Not currently used.
--------------------------------------------------------------------------------
Returns:
Returns tuple (k,eigval):
* k:        Not currently used
* eigval:   Spectral radius
--------------------------------------------------------------------------------
"""
function spectralR(hj,X,P,rho)
  k = true
  J = hj(X,P)
  Y = Preconditioner(hj,X,P,jac="User")
  YJ = abs.(Y)*diam.(J)/2
  if (length(YJ)>1)
    eigval = eigmax(YJ)
  else
    eigval = YJ[1]
  end
  return k,eigval
end
=#
#=
"""
--------------------------------------------------------------------------------
Function: BoundaryTest
--------------------------------------------------------------------------------
Description:
Runs novel test for inclusion of implicit function
--------------------------------------------------------------------------------
Inputs:
* h:                function - h(x,p) potential implicit function
* hj:               function - Jacobian wrt x of potential implicit function
* X0:               Vector{Interval{Float64}} - Bounds for dependent variables
* P:                Vector{Interval{Float64}} - Bounds for independent variables
* opt:              Array{Any,1} - Parameters used
* PIcert:           Bool: Partial certification flag
* PIflag:           Bool: Partial inclusion flag
* Iflag:            Bool: Inclusion flag
* Eflag:            Bool: Exclusion flag
* inclusionLow:     Vector{Bool}: Array of boolean for inclusion test
* inclusionHigh:    Vector{Bool}:  Array of boolean for inclusion test
--------------------------------------------------------------------------------
Returns:
Returns tuple (PIcert,PIflag,Iflag,Eflag):
* PIcert:           Bool: Partial certification flag
* PIflag:           Bool: Partial inclusion flag
* Iflag:            Bool: Inclusion flag
* Eflag:            Bool: Exclusion flag
--------------------------------------------------------------------------------
"""
function BoundaryTest(h,hj,X0,X,P,opt,PIcert,PIflag,Iflag,Eflag,inclusionLow,inclusionHigh)
  PIflagT = PIflag
  PIflagB = PIflag
  Xtemp = copy(X)
  Xin = copy(X)
  Pin = copy(P)
  ku,PIflagT = partialIncTop(h,Xin,Pin,PIflagT,inclusionHigh)
  kl,PIflagB = partialIncBot(h,Xin,Pin,PIflagB,inclusionLow)
  rho = 1.01
  eDflag = false
  if (Strict_XinY(X,X0) || ((ku==-1) && (kl==-1)))
    blank,rho = spectralR(hj,Xin,Pin,rho)
    PIcert = true
    PIflag = false
    if (rho<1.0)
      Pmid = mid.(Pin)
      if (opt[4]=="NewtonGS")
        Xnew,Xtemp,Eflag,Iflag,eDflag = PI_NewtonGS(Xin,Pmid,hj,h,opt,Eflag,Iflag,eDflag)
      elseif (opt[4]=="KrawczykCW")
        Xnew,Eflag,Iflag,blank1,blank2 =  MC_KrawczykCW(Xin,Pmid,hj,h,opt,Eflag,Iflag)
      end
    else
      Iflag = false
    end
  end
  return PIcert,PIflag,Iflag,Eflag
end
=#
#=
function NewBoundaryTest(Xtemp,X,P,Iflag)
  nx = length(X)
  kmax = 1000
  k = 1
  l = 0
  ztemp = 0.0
  eps = 1.0E-12
  Xf = copy(X)
  Xnew = copy(X)
  Pf = copy(P)
  l = SingleStep(Xnew,Pf)
  if (l == 0)
    k = kmax
  end
  for i = 1:nx
    r[2*i-1] = Xf[i].lo - Xnew[i].lo + eps
    r[2*i] = Xnew[i].hi - Xf[i].hi + eps
    z[2*i-1] = Xf[i].lo
    z[2*i] = Xf[i].hi
  end
  # construct Jacobian
  Jr =

end

function SingleStep()
end
=#

ext: jl


text:
function convert_wildcards!(account, parameters)
    for r = 1:3
        while account.wc_count[r] >= parameters.wc_upgrade_rate[r] && account.wc_count[r] >= parameters.wc_upgrade_threshold[r]
            account.wc_count[r] -= parameters.wc_upgrade_rate[r]
            account.wc_count[r+1] += 1
        end
    end
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ eb97df10-76f1-11eb-21da-378e76583695
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)

	using MLJ
	using CSV
	using DataFrames
	using PlutoUI
	using Test
	using Printf
	# using Plots
end

# ╔═╡ 299e3eaa-7709-11eb-2b3e-13dcfbd9f450
begin
	include("./utils.jl");
	include("./text_proc.jl");
end

# ╔═╡ 9cf87de4-76f1-11eb-221a-8fbc8dad79fd
md"""
## C03/w01: Predicting sentiment from product reviews

The goal of this first notebook is to explore logistic regression and feature engineering with Julia and MLJ.

In this notebook we will use product review data from Amazon.com to predict whether the sentiments about a product (from its reviews) are positive or negative.

  -  Use DataFrames to do some feature engineering
  - Train a logistic regression model to predict the sentiment of product reviews.
  -  Inspect the weights (coefficients) of a trained logistic regression model.
  - Make a prediction (both class and probability) of sentiment for a new product review.
  -  Given the logistic regression weights, predictors and ground truth labels, write a function to compute the **accuracy** of the model.
  -  Inspect the coefficients of the logistic regression model and interpret their meanings.
  -  Compare multiple logistic regression models.
"""

# ╔═╡ fea86e6c-76f1-11eb-0188-a1d8b7bcf47b
md"""
### Data Preparation

We will use a dataset consisting of baby product reviews on Amazon.com. Let us load it and see a preview of what it looks like.
"""

# ╔═╡ 17de9582-76f2-11eb-03cf-43ad5b07d120
begin
	products = train = CSV.File("../../ML_UW_Spec/C03/data/amazon_baby.csv";
		header=true) |> DataFrame;

	size(products)
end

# ╔═╡ f29e95dc-7d50-11eb-1ecb-51f578ef29f0
md"""
###### Deal with missing values

  1. Filling missing review with empty string
  2. drop other missing rows
"""

# ╔═╡ 01417a14-7d51-11eb-394e-a36a82d4a40a
replace!(products.review, missing =>"");

# ╔═╡ ade2915c-7709-11eb-2fce-45999f1e9ffa
begin
	## drop missing value if any...
	before_rm_na = size(products)
	dropmissing!(products);
	before_rm_na, size(products)
end

# ╔═╡ 67034e64-76f2-11eb-2a25-c53adf9adde6
md"""
### Build the word count vector for each review

Let us explore a specific example of a baby product.
"""

# ╔═╡ 7680d44c-76f2-11eb-1ffe-f92fdd5e42b7
products[270, :]

# ╔═╡ db3bbb6e-76f6-11eb-09c1-e7f0392c634d
products[70, :review]

# ╔═╡ 76610114-76f2-11eb-3e22-71a763599e47
md"""
Now, we will perform 2 simple data transformations:

  - Remove punctuation (using a regular expression) and
  - Transform the reviews into word-counts.

Aside. In this notebook, we remove all punctuations for the sake of simplicity. A smarter approach to punctuations would preserve phrases such as "I'd", "would've", "hadn't" and so forth.

"""

# ╔═╡ 762adc54-76f2-11eb-259d-8f6cae34ef20
remove_punctuation(" Voici un example; ponctuation: elimine? Vraiment?? Oui! *****")

# ╔═╡ 75e24a66-76f2-11eb-2711-f53a6c04fb99
products[!, :review] = remove_punctuation.(products[:, :review]);

# ╔═╡ 75bfd17a-76f2-11eb-3285-31b3eaa96a62
products[70, :review]

# ╔═╡ 3b78e5bc-7707-11eb-1e32-2116b3156727
remove_punctuation("   Voici un example; ponctuation: elimine. Vraiment?? Oui! *****")   |> word_count

# ╔═╡ b398fef4-7706-11eb-2b36-17f85f356a3c
first(products, 3)

# ╔═╡ 42d86fd8-7878-11eb-0e98-c73dd498e285
md"""
#### Clean up review
"""

# ╔═╡ 861170cc-7877-11eb-00c7-0f50fe4ceeea
## clean review
begin
	products[!, :review_clean] = cleanup.(products.review);
	select!(products, Not(:review));
	size(products)
end

# ╔═╡ d4844b50-770b-11eb-221c-0f281a8d8eec
md"""
#### Extract sentiment

We will ignore all reviews with `rating=3`, since they tend to have a neutral sentiment.

Now, we will assign reviews with a rating of 4 or higher to be positive reviews, while the ones with rating of 2 or lower are negative.
For the sentiment column, we use +1 for the positive class label and -1 for the negative class label.

"""

# ╔═╡ f40d002a-770b-11eb-3d91-8916193cf7a3
begin
	filter!(r -> r.rating != 3, products);
	products.sentiment = ifelse.(products.rating .> 3, 1 ,-1);
	size(products)
end

# ╔═╡ ce050256-7aee-11eb-1474-b96dcad21315
md"""
### MLJ Data Preparation
"""

# ╔═╡ 5fcd5f68-7715-11eb-3667-7d4960547cff
md"""
#### CountVectorizer
"""

# ╔═╡ 187d5a38-7875-11eb-01c1-bbc44e1d15c0
token_ix = gen_vocab.(products.review_clean) |>
   vtext -> reduce(vcat, vtext) |>
   unique |>
   words -> filter(w -> length(w) > 2, words) |>
   vocab -> Dict{String, Int}(w => ix for (ix, w) in enumerate(vocab))
## 58836 unique words -> 58584

# ╔═╡ 0f7839a6-787e-11eb-2887-f34311d11b56
length(token_ix)

# ╔═╡ a6d8c4c2-771c-11eb-30bc-332df04ade43
transform!(products, :review_clean => h_encode => :word_count);

# ╔═╡ 4dda42ea-774b-11eb-2136-b345e017b297
names(products)

# ╔═╡ 9c2dd014-7751-11eb-1ecb-51f578ef29f0
first(products, 3)

# ╔═╡ e265712c-7aee-11eb-28c8-9952ca7fbc46
## before coercion
schema(products)

# ╔═╡ 375e9688-7af2-11eb-08fe-efb6d6b91863
begin
	## Note we only use :word_count as a feature for our Logistic Regr. classifier
	## ST = Scientific Type

	const Features_ST = Dict{Symbol, DataType}(
		:word_count => Continuous,             ## num. of unique words...
	# 	:word_count => Count,
	)

	Target_ST = Dict{Symbol, DataType}(
		:sentiment => Multiclass{2}
	)
end

# ╔═╡ ef617988-7d59-11eb-2f81-69e21d6c7ffa
elscitype(products.word_count)

# ╔═╡ 69aa0492-7d5c-11eb-2f81-69e21d6c7ffa
scitype(products.word_count)

# ╔═╡ 27f66dec-7d4f-11eb-29a5-d79c586de1d8
## Coerce
begin
	coerce_map = [Features_ST..., Target_ST...]
	coerce!(products, coerce_map...);

	# then
	# coerce!(products, autotype(products, :discrete_to_continuous))
	schema(products)
end

# ╔═╡ eced8604-7d44-11eb-246e-a77066a91eb1
# reduce(vcat, products[!, :name]) |> unique |> length # == 30628
# name => Multiclass{N} ?

# ╔═╡ 7e24bb38-770a-11eb-03ef-590069fda67c
md"""
#### Split data into train and test sets
"""

# ╔═╡ b3580318-7706-11eb-0b23-6b817d62e154
begin
	train_data, test_data = train_test_split(products; split=0.8, seed=42)
	size(train_data), size(test_data)
end

# ╔═╡ 909e2708-770b-11eb-0955-6f3c7c6d19ec
names(train_data)

# ╔═╡ b33d33f8-7706-11eb-1b83-130c4fb81eb2
md"""
### Train a sentiment classifier with logistic regression

We will now use logistic regression to create a sentiment classifier on the training data. This model will use the column word_count as a feature and the column sentiment as the target.
"""

# ╔═╡ a572e232-770a-11eb-12ce-d5cf2e70d878
@load LogisticClassifier pkg=MLJLinearModels

# ╔═╡ e3809494-7af6-11eb-33b9-4d503e9aa1fc
function data_prep_(df::DF; features=Features_ST,
	target=Target_ST) where {DF <: Union{DataFrame, SubDataFrame}}

 	y_target = collect(keys(target))[1]
    x_features = collect(keys(features))

	X = select(df, x_features)
	y = select(df, y_target)
	y = y[:, y_target]

   	return (X, y)
end

# ╔═╡ 29716c86-7d53-11eb-1de9-7ba4e5be7112
begin
	Xₜ, yₜ = data_prep_(train_data)
	schema(Xₜ)
end

# ╔═╡ be5a7876-7d56-11eb-145f-45c61470f7cc
first(Xₜ, 3)

# ╔═╡ 46c01a6c-7d53-11eb-0e1b-cbc471d0a34e
typeof(Xₜ), typeof(yₜ)

# ╔═╡ a556fdf6-770a-11eb-18c3-3fc37fe4fe32
function run_logistic_regression(X, y)
	mach = machine(MLJLinearModels.LogisticClassifier(), X, y)
	fit!(mach) # , rows=1:nrows, verbosity=1, force=false)
	mach
end

# ╔═╡ b30dea8a-7706-11eb-266b-eb9027195a0b
begin
	sentiment_mdl = run_logistic_regression(Xₜ, yₜ);
	fitted_params(sentiment_mdl) # report(sentiment_mdl)
	# ref: https://github.com/alan-turing-institute/MLJ.jl/issues/492
end

# ╔═╡ fc20085e-7d53-11eb-2f81-69e21d6c7ffa

# ╔═╡ Cell order:
# ╟─9cf87de4-76f1-11eb-221a-8fbc8dad79fd
# ╠═eb97df10-76f1-11eb-21da-378e76583695
# ╟─fea86e6c-76f1-11eb-0188-a1d8b7bcf47b
# ╠═17de9582-76f2-11eb-03cf-43ad5b07d120
# ╟─f29e95dc-7d50-11eb-1ecb-51f578ef29f0
# ╠═01417a14-7d51-11eb-394e-a36a82d4a40a
# ╠═ade2915c-7709-11eb-2fce-45999f1e9ffa
# ╠═299e3eaa-7709-11eb-2b3e-13dcfbd9f450
# ╟─67034e64-76f2-11eb-2a25-c53adf9adde6
# ╠═7680d44c-76f2-11eb-1ffe-f92fdd5e42b7
# ╠═db3bbb6e-76f6-11eb-09c1-e7f0392c634d
# ╟─76610114-76f2-11eb-3e22-71a763599e47
# ╠═762adc54-76f2-11eb-259d-8f6cae34ef20
# ╠═75e24a66-76f2-11eb-2711-f53a6c04fb99
# ╠═75bfd17a-76f2-11eb-3285-31b3eaa96a62
# ╠═3b78e5bc-7707-11eb-1e32-2116b3156727
# ╠═b398fef4-7706-11eb-2b36-17f85f356a3c
# ╟─42d86fd8-7878-11eb-0e98-c73dd498e285
# ╠═861170cc-7877-11eb-00c7-0f50fe4ceeea
# ╟─d4844b50-770b-11eb-221c-0f281a8d8eec
# ╠═f40d002a-770b-11eb-3d91-8916193cf7a3
# ╟─ce050256-7aee-11eb-1474-b96dcad21315
# ╟─5fcd5f68-7715-11eb-3667-7d4960547cff
# ╠═187d5a38-7875-11eb-01c1-bbc44e1d15c0
# ╠═0f7839a6-787e-11eb-2887-f34311d11b56
# ╠═a6d8c4c2-771c-11eb-30bc-332df04ade43
# ╠═4dda42ea-774b-11eb-2136-b345e017b297
# ╠═9c2dd014-7751-11eb-1ecb-51f578ef29f0
# ╠═e265712c-7aee-11eb-28c8-9952ca7fbc46
# ╠═375e9688-7af2-11eb-08fe-efb6d6b91863
# ╠═ef617988-7d59-11eb-2f81-69e21d6c7ffa
# ╠═69aa0492-7d5c-11eb-2f81-69e21d6c7ffa
# ╠═27f66dec-7d4f-11eb-29a5-d79c586de1d8
# ╠═eced8604-7d44-11eb-246e-a77066a91eb1
# ╟─7e24bb38-770a-11eb-03ef-590069fda67c
# ╠═b3580318-7706-11eb-0b23-6b817d62e154
# ╠═909e2708-770b-11eb-0955-6f3c7c6d19ec
# ╟─b33d33f8-7706-11eb-1b83-130c4fb81eb2
# ╠═a572e232-770a-11eb-12ce-d5cf2e70d878
# ╠═e3809494-7af6-11eb-33b9-4d503e9aa1fc
# ╠═29716c86-7d53-11eb-1de9-7ba4e5be7112
# ╠═be5a7876-7d56-11eb-145f-45c61470f7cc
# ╠═46c01a6c-7d53-11eb-0e1b-cbc471d0a34e
# ╠═a556fdf6-770a-11eb-18c3-3fc37fe4fe32
# ╠═b30dea8a-7706-11eb-266b-eb9027195a0b
# ╠═fc20085e-7d53-11eb-2f81-69e21d6c7ffa


text:
using Pkg, Conda

# this line doesn't seem to be hit
println("building ECCOtour package")

#Pkg.add(url="https://github.com/ggebbie/SigmaShift.jl")
#Pkg.instantiate()

if lowercase(get(ENV, "CI", "false")) == "true"

    ENV["PYTHON"] = ""
    Pkg.build("PyCall")

    Conda.add("matplotlib",channel="conda-forge")
    Conda.add("shapely",channel="conda-forge")
    Conda.add("cartopy=0.20.0",channel="conda-forge")

end


text:
"""
    append_axes(x, y)

Returns the axes for `append(x, y)`.
"""
append_axes(x::AbstractVector, y::AbstractVector) = (append_axis(axes(x, 1), axes(y, 1)),)

"""
    append_axis(x, y)

Returns the appended axes `x` and `y`. New subtypes of `AbstractIndex` must
implement a unique `append_axis` method.

## Examples
```jldoctest
julia> using AbstractIndices: append_axis!

julia> x, y = Index(UnitMRange(1, 10)), SimpleIndex(UnitMRange(1, 10));

julia> append_axis(x, y)
Index(UnitMRange(1:20) => OneToMRange(20))

julia> append_axis(y, x)
SimpleIndex(UnitMRange(1:20))
```
"""
function append_axis(x::Index, y::Index)
    return Index{combine_names(x, y)}(append_keys(x, y), append_values(x, y))
end
function append_axis(x::SimpleIndex, y::SimpleIndex)
    return SimpleIndex{append_names(x, y)}(append_values(x, y))
end
append_axis(x::AbstractIndex, y::AbstractVector) = append_axis(promote(x, y)...)
append_axis(x::AbstractVector, y::AbstractIndex) = append_axis(promote(x, y)...)

"""
    append_keys(x, y)

Returns the appropriate keys of and index within the operation `append_axis(x, y)`

See also: [`append_axis`](@ref)
"""
append_keys(x, y) = cat_keys(x, y)

"""
    append_values(x, y)

Returns the appropriate values of and index within the operation `append_axis(x, y)`

See also: [`append_axis`](@ref)
"""
append_values(x, y) = cat_values(x, y)

"""
    append_axes!(x, y)

Returns the axes for `append!(x, y)`.
"""
append_axes!(x::AbstractVector, y::AbstractVector) = (append_axis!(axes(x, 1), axes(y, 1)),)

"""
    append_axis!(x, y)

Returns the appended axes `x` and `y`. New subtypes of `AbstractIndex` must
implement a unique `append_axis!` method.

## Examples
```jldoctest
julia> using AbstractIndices: append_axis!

julia> x, y = Index(UnitMRange(1, 10)), SimpleIndex(UnitMRange(1, 10));

julia> append_axis!(x, y)
Index(UnitMRange(1:20) => OneToMRange(20))

julia> append_axis!(y, x)
SimpleIndex(UnitMRange(1:30))
```
"""
function append_axis!(x::Index, y)
    _append_keys!(keys(x), y)
    set_length!(values(x), length(x) + length(y))
    return x
end
function append_axis!(x::SimpleIndex, y)
    set_length!(x, length(x) + length(y))
    return x
end

_append_keys!(x, y) = __append_keys!(Continuity(x), x, y)
__append_keys!(::ContinuousTrait, x, y) = set_length!(x, length(x) + length(y))
__append_keys!(::DiscreteTrait, x, y) = make_unique!(x, keys(y))


text:
@testset "Rewrite Rules" begin
    a, b = LogicalSymbol.([:a, :b])

end

@testset "Propositional Proofs" begin
    a, b, c, d = LogicalSymbol.([:a, :b, :c, :d])

    # trivially true, since we assert any proposition as true when passed as a parameter of `tableau`
    @test tableau(a)
    @test tableau(b)
    @test tableau(¬a)
    @test tableau(a ∧ b)
    @test tableau(a ∨ b)
    @test tableau(a → b)
    @test tableau(a ⟷ b)

    # multiple premises, still trivially true
    @test tableau([a, b])
    @test tableau(a, b)  # alternate preferable syntax
    @test tableau(a → b, a, b)

    # some contradictory premises should yield false
    @test !tableau(a, ¬a)
    @test !tableau(a ∧ b, ¬a)
    @test !tableau(a → b, a, ¬b)
    @test !tableau(a ⟷ b, ¬a, b)
    @test !tableau(a ∧ b ∨ c, ¬c, ¬a)
    @test !tableau(a → b, b → c, c → d, a, ¬d)
end

@testset "Assertion Proofs" begin
    @testset "Rule Combinations" begin
        @symbols a b c d
        @symbols p q

        conjunction = rule_by_name(PropositionalCalculus, "Conjunction")
        premises1 = Set{AbstractExpression}([a, b])
        premises2 = Set{AbstractExpression}([a, b, c])
        premises3 = Set{AbstractExpression}([a, b, c, d])

        @test rule_combinations(conjunction, premises1) == Set([
            Deductive.SymbolMap(p => a, q => a),
            Deductive.SymbolMap(p => a, q => b),
            Deductive.SymbolMap(p => b, q => a),
            Deductive.SymbolMap(p => b, q => b),
        ])
        @test rule_combinations(conjunction, premises2) == Set([
            Deductive.SymbolMap(p => a, q => a),
            Deductive.SymbolMap(p => a, q => b),
            Deductive.SymbolMap(p => a, q => c),
            Deductive.SymbolMap(p => b, q => a),
            Deductive.SymbolMap(p => b, q => b),
            Deductive.SymbolMap(p => b, q => c),
            Deductive.SymbolMap(p => c, q => a),
            Deductive.SymbolMap(p => c, q => b),
            Deductive.SymbolMap(p => c, q => c)
        ])
        @test length(rule_combinations(conjunction, premises3)) == 16
        #  lol did u really think I'd write out all 16 combinations 🤣🤣🤣

        modus_ponens = rule_by_name(PropositionalCalculus, "Modus Ponens")
        premises4 = Set{AbstractExpression}([
            a → b,
            a
        ])
        premises5 = Set{AbstractExpression}([
            a → b,
            (a → b) → c
        ])

        @test rule_combinations(modus_ponens, premises4) == Set([
            Deductive.SymbolMap(p => a, q => b)
        ])
        @test rule_combinations(modus_ponens, premises5) == Set([
            Deductive.SymbolMap(p => a → b, q => c)
        ])

        double_negation = rule_by_name(PropositionalCalculus, "Double Negation Introduction")
        premises6 = Set{AbstractExpression}([
            a,
            b,
            a ∧ b
        ])

        @test rule_combinations(double_negation, premises6) == Set([
            Deductive.SymbolMap(p => a),
            Deductive.SymbolMap(p => b),
            Deductive.SymbolMap(p => a ∧ b)
        ])
    end
end

# @testset "Predicate Proofs" begin
#     a, b = LogicalSymbol.([:a, :b])
#     x, y, z = FreeVariable.([:x, :y, :z])
#     P, Q = Predicate.([:P, :Q])

#     # test cases taken from https://www.esf.kfi.zcu.cz/logika/opory/ia008/tableau_pred-sol.pdf

#     ex_a = Ā(x, P(x)) ⟷ ¬Ē(x, ¬P(x))
#     @test !tableau(¬ex_a)

#     ex_b = Ā(x, P(x) → Q(x)) → (Ā(x, P(x)) → Ā(x, Q(x)))
#     @test !tableau(¬ex_b)

#     ex_c = Ā(x, P(x) ∧ Q(x)) ⟷ (Ā(x, P(x)) ∧ Ā(x, Q(x)))
#     @test !tableau(¬ex_c)

#     # support for multiple free variables later...
#     # ex_d = Ē(y, Ā(x, P(x, y) ⟷ P(x, x))) → ¬Ā(x, Ē(y, Ā(z, P(z, y) ⟷ P(z, x))))
#     # @test !tableau(¬ex_d)
# end


text:
module ModCartesianIndicesBase
export ModCartesianIndices

using Base: tail
import Base: getindex, length, size, IndexStyle, IndexCartesian, in, iterate, mod,
    first, last
@inline mod_tuple(c::NTuple{N,Int}, size::NTuple{N,Int}) where N  =
    ntuple(k->mod(c[k]-1,size[k])+1, Val(N))

@inline mod(c::CartesianIndex{N}, size::NTuple{N,Int}) where N =
    CartesianIndex(mod_tuple(c.I, size))

@inline add_offset_mod(i::NTuple{N,Int}, start::NTuple{N,Int}, size::NTuple{N,Int}, periodic::NTuple{N,Bool}) where N =
    ntuple(k->( periodic[k]  ? add_offset_mod(i[k], start[k], size[k]) : i[k] ), Val(N))

@inline add_offset_mod(i::Int, start::Int, size::Int) =
    start + mod(i-start, size)

struct ModCartesianIndices{N,K} <:AbstractArray{CartesianIndex{N},N}
    size::NTuple{N,Int}
    iter::CartesianIndices{N,K}
    start::CartesianIndex{N}
    stop::CartesianIndex{N}
    periodic::NTuple{N,Bool}
end

getindex(A::ModCartesianIndices{N}, I::Vararg{Int,N}) where N = mod(getindex(A.iter,I...),A.size)
IndexStyle(A::ModCartesianIndices) = IndexCartesian()
function ModCartesianIndices(size::NTuple{N,Int}, iranges::UnitRange...) where {N}
    start = CartesianIndex(ntuple(k->first(iranges[k]), Val(N)))
    stop  = CartesianIndex(ntuple(k->last(iranges[k]), Val(N)))
    ModCartesianIndices(size, CartesianIndices(iranges), start, stop, ntuple(k->true,Val(N)))
end

function ModCartesianIndices(size::NTuple{N,Int}, start::CartesianIndex{N}, stop::CartesianIndex,  periodic::NTuple{N,Bool}=ntuple(k->true,Val(N))) where {N}
    iranges = ntuple(k-> ((start[k] <= stop[k]) ? (start[k]:stop[k]) : (start[k]:(stop[k]+size[k]))), Val(N))
    ModCartesianIndices(size, CartesianIndices(iranges), start, stop, periodic)
end

@inline size(m::ModCartesianIndices) = size(m.iter)
@inline length(m::ModCartesianIndices) = length(m.iter)

struct ModUnitRange{I} <:AbstractUnitRange{I}
    size::I
    iter::UnitRange{I}
    periodic::Bool

    function ModUnitRange(size::I, iter::UnitRange{I}, periodic=true) where I<:Integer
        if periodic
            e = last(iter)
            l = length(iter)
            e_aligned = mod(e-1,size)+1
            f_aligned = e_aligned - l + 1
            new{I}(size, f_aligned:e_aligned, periodic)
        else
            new{I}(size, max(1,start(iter)):min(size,last(iter)), periodic)
        end
    end
end

first(r::ModUnitRange) = mod(first(r.iter)-1,r.size)+1
last(r::ModUnitRange) = mod(last(r.iter)-1,r.size)+1

function in(i::Integer, r::ModUnitRange{<:Integer})
    if r.periodic
        i_aligned =  mod(i-1,r.size)+1
        return (i_aligned ∈ r.iter) | (i_aligned ∈ r.iter .+ r.size)
    else
        return in(i, r.iter)
    end
end

function iterate(m::ModUnitRange)
    tuple = iterate(m.iter)
    if tuple != nothing
        index, state = tuple
        mod(index-1, m.size)+1, state
    end
end

function iterate(m::ModUnitRange, state)
    tuple = iterate(m.iter, state)
    if tuple != nothing
        index, state = tuple
        mod(index-1, m.size)+1, state
    end
end

@inline length(m::ModUnitRange) = length(m.iter)

nbindexlist(index::Base.CartesianIndex{N}, size::NTuple{N,Int}, periodic=ntuple(k->false, Val(N))) where {N}  =
    ModCartesianIndices(size, index+CartesianIndex(ntuple(k->-1, Val(N))), index+CartesianIndex(ntuple(k->1, Val(N))), periodic)
nbindexlist(index::Int, size::NTuple{1,Int}, periodic=false)  =
    ModUnitRange(size[1], index-1:index+1, periodic)

end # module


text:
int_rules_1_1_3_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.3.6*(g*x)^m*(a+b*x^n)^p*(c+d*x^n)^q*(e+f*x^n)^r =#
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => ((~g) ^ ~m / (~n * (~b) ^ (Simplify((~m + 1) / ~n) - 1))) * Subst(Antiderivative((~b * ~x) ^ ((~p + Simplify((~m + 1) / ~n)) - 1) * (~c + ~d * ~x) ^ ~q * (~e + ~f * ~x) ^ ~r, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && ((IntegerQ(~m) || GtQ(~g, 0)) && IntegerQ(Simplify((~m + 1) / ~n)))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~(b') * (~x) ^ ~(n')) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (((~g) ^ ~m * (~b) ^ IntPart(~p) * (~b * (~x) ^ ~n) ^ FracPart(~p)) / (~x) ^ (~n * FracPart(~p))) * Antiderivative((~x) ^ (~m + ~n * ~p) * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && ((IntegerQ(~m) || GtQ(~g, 0)) && Not(IntegerQ(Simplify((~m + 1) / ~n))))
    @apply_utils Antiderivative((~g * ~x) ^ ~m * (~(b') * (~x) ^ ~(n')) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (((~g) ^ IntPart(~m) * (~g * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => Antiderivative(ExpandIntegrand((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n], ~x) && (IGtQ(~p, -2) && (IGtQ(~q, 0) && IGtQ(~r, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (1 / ~n) * Subst(Antiderivative((~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q * (~e + ~f * ~x) ^ ~r, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~q, ~r], ~x) && EqQ((~m - ~n) + 1, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => Antiderivative((~x) ^ (~m + ~n * (~p + ~q + ~r)) * (~b + ~a * (~x) ^ -(~n)) ^ ~p * (~d + ~c * (~x) ^ -(~n)) ^ ~q * (~f + ~e * (~x) ^ -(~n)) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (IntegersQ(~p, ~q, ~r) && NegQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q * (~e + ~f * ~x) ^ ~r, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~q, ~r], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~g * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (((~g) ^ IntPart(~m) * (~g * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => With([k = GCD(~m + 1, ~n)], (1 / k) * Subst(Antiderivative((~x) ^ ((~m + 1) / k - 1) * (~a + ~b * (~x) ^ (~n / k)) ^ ~p * (~c + ~d * (~x) ^ (~n / k)) ^ ~q * (~e + ~f * (~x) ^ (~n / k)) ^ ~r, ~x), ~x, (~x) ^ k) <-- k != 1) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r], ~x) && (IGtQ(~n, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~q * (~e + ~(f') * (~x) ^ ~n) ^ ~r, ~x) => With([k = Denominator(~m)], (k / ~g) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + (~b * (~x) ^ (k * ~n)) / (~g) ^ ~n) ^ ~p * (~c + (~d * (~x) ^ (k * ~n)) / (~g) ^ ~n) ^ ~q * (~e + (~f * (~x) ^ (k * ~n)) / (~g) ^ ~n) ^ ~r, ~x), ~x, (~g * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p, ~q, ~r], ~x) && (IGtQ(~n, 0) && FractionQ(~m))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (-((~b * ~e - ~a * ~f)) * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q) / (~a * ~b * ~g * ~n * (~p + 1)) + (1 / (~a * ~b * ~n * (~p + 1))) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q - 1) * Simp(~c * (~b * ~e * ~n * (~p + 1) + (~b * ~e - ~a * ~f) * (~m + 1)) + ~d * (~b * ~e * ~n * (~p + 1) + (~b * ~e - ~a * ~f) * (~m + ~n * ~q + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m], ~x) && (IGtQ(~n, 0) && (LtQ(~p, -1) && (GtQ(~q, 0) && Not(EqQ(~q, 1) && SimplerQ(~b * ~c - ~a * ~d, ~b * ~e - ~a * ~f)))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~q * (~e + ~(f') * (~x) ^ ~n), ~x) => ((~g) ^ (~n - 1) * (~b * ~e - ~a * ~f) * (~g * ~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q + 1)) / (~b * ~n * (~b * ~c - ~a * ~d) * (~p + 1)) - ((~g) ^ ~n / (~b * ~n * (~b * ~c - ~a * ~d) * (~p + 1))) * Antiderivative((~g * ~x) ^ (~m - ~n) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q * Simp(~c * (~b * ~e - ~a * ~f) * ((~m - ~n) + 1) + (~d * (~b * ~e - ~a * ~f) * (~m + ~n * ~q + 1) - ~b * ~n * (~c * ~f - ~d * ~e) * (~p + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~q], ~x) && (IGtQ(~n, 0) && (LtQ(~p, -1) && GtQ((~m - ~n) + 1, 0)))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~q * (~e + ~(f') * (~x) ^ ~n), ~x) => (-((~b * ~e - ~a * ~f)) * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q + 1)) / (~a * ~g * ~n * (~b * ~c - ~a * ~d) * (~p + 1)) + (1 / (~a * ~n * (~b * ~c - ~a * ~d) * (~p + 1))) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q * Simp(~c * (~b * ~e - ~a * ~f) * (~m + 1) + ~e * ~n * (~b * ~c - ~a * ~d) * (~p + 1) + ~d * (~b * ~e - ~a * ~f) * (~m + ~n * (~p + ~q + 2) + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~q], ~x) && (IGtQ(~n, 0) && LtQ(~p, -1))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (~e * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q) / (~a * ~g * (~m + 1)) - (1 / (~a * (~g) ^ ~n * (~m + 1))) * Antiderivative((~g * ~x) ^ (~m + ~n) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ (~q - 1) * Simp(~c * (~b * ~e - ~a * ~f) * (~m + 1) + ~e * ~n * (~b * ~c * (~p + 1) + ~a * ~d * ~q) + ~d * ((~b * ~e - ~a * ~f) * (~m + 1) + ~b * ~e * ~n * (~p + ~q + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (IGtQ(~n, 0) && (GtQ(~q, 0) && (LtQ(~m, -1) && Not(EqQ(~q, 1) && SimplerQ(~e + ~f * (~x) ^ ~n, ~c + ~d * (~x) ^ ~n)))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (~f * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q) / (~b * ~g * (~m + ~n * (~p + ~q + 1) + 1)) + (1 / (~b * (~m + ~n * (~p + ~q + 1) + 1))) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ (~q - 1) * Simp(~c * ((~b * ~e - ~a * ~f) * (~m + 1) + ~b * ~e * ~n * (~p + ~q + 1)) + (~d * (~b * ~e - ~a * ~f) * (~m + 1) + ~f * ~n * ~q * (~b * ~c - ~a * ~d) + ~b * ~e * ~d * ~n * (~p + ~q + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p], ~x) && (IGtQ(~n, 0) && (GtQ(~q, 0) && Not(EqQ(~q, 1) && SimplerQ(~e + ~f * (~x) ^ ~n, ~c + ~d * (~x) ^ ~n))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (~f * (~g) ^ (~n - 1) * (~g * ~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q + 1)) / (~b * ~d * (~m + ~n * (~p + ~q + 1) + 1)) - ((~g) ^ ~n / (~b * ~d * (~m + ~n * (~p + ~q + 1) + 1))) * Antiderivative((~g * ~x) ^ (~m - ~n) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * Simp(~a * ~f * ~c * ((~m - ~n) + 1) + (~a * ~f * ~d * (~m + ~n * ~q + 1) + ~b * (~f * ~c * (~m + ~n * ~p + 1) - ~e * ~d * (~m + ~n * (~p + ~q + 1) + 1))) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p, ~q], ~x) && (IGtQ(~n, 0) && GtQ(~m, ~n - 1))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (~e * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q + 1)) / (~a * ~c * ~g * (~m + 1)) + (1 / (~a * ~c * (~g) ^ ~n * (~m + 1))) * Antiderivative((~g * ~x) ^ (~m + ~n) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * Simp(((~a * ~f * ~c * (~m + 1) - ~e * (~b * ~c + ~a * ~d) * (~m + ~n + 1)) - ~e * ~n * (~b * ~c * ~p + ~a * ~d * ~q)) - ~b * ~e * ~d * (~m + ~n * (~p + ~q + 2) + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p, ~q], ~x) && (IGtQ(~n, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative(((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~e + ~(f') * (~x) ^ ~n)) / (~c + ~(d') * (~x) ^ ~n), ~x) => Antiderivative(ExpandIntegrand(((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~e + ~f * (~x) ^ ~n)) / (~c + ~d * (~x) ^ ~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => ~e * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) + (~f / (~e) ^ ~n) * Antiderivative((~g * ~x) ^ (~m + ~n) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p, ~q], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => ~e * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ (~r - 1), ~x) + (~f / (~e) ^ ~n) * Antiderivative((~g * ~x) ^ (~m + ~n) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ (~r - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p, ~q], ~x) && (IGtQ(~n, 0) && IGtQ(~r, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => -(Subst(Antiderivative(((~a + ~b * (~x) ^ -(~n)) ^ ~p * (~c + ~d * (~x) ^ -(~n)) ^ ~q * (~e + ~f * (~x) ^ -(~n)) ^ ~r) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r], ~x) && (ILtQ(~n, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => With([k = Denominator(~m)], (-k / ~g) * Subst(Antiderivative(((~a + ~b * (~g) ^ -(~n) * (~x) ^ (-k * ~n)) ^ ~p * (~c + ~d * (~g) ^ -(~n) * (~x) ^ (-k * ~n)) ^ ~q * (~e + ~f * (~g) ^ -(~n) * (~x) ^ (-k * ~n)) ^ ~r) / (~x) ^ (k * (~m + 1) + 1), ~x), ~x, 1 / (~g * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p, ~q, ~r], ~x) && (ILtQ(~n, 0) && FractionQ(~m))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => -((~g * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative(((~a + ~b * (~x) ^ -(~n)) ^ ~p * (~c + ~d * (~x) ^ -(~n)) ^ ~q * (~e + ~f * (~x) ^ -(~n)) ^ ~r) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p, ~q, ~r], ~x) && (ILtQ(~n, 0) && Not(RationalQ(~m)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * (~x) ^ (k * ~n)) ^ ~p * (~c + ~d * (~x) ^ (k * ~n)) ^ ~q * (~e + ~f * (~x) ^ (k * ~n)) ^ ~r, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p, ~q, ~r], ~x) && FractionQ(~n)
    @apply_utils Antiderivative((~g * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (((~g) ^ IntPart(~m) * (~g * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p, ~q, ~r], ~x) && FractionQ(~n)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ Simplify(~n / (~m + 1))) ^ ~p * (~c + ~d * (~x) ^ Simplify(~n / (~m + 1))) ^ ~q * (~e + ~f * (~x) ^ Simplify(~n / (~m + 1))) ^ ~r, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~q, ~r], ~x) && IntegerQ(Simplify(~n / (~m + 1)))
    @apply_utils Antiderivative((~g * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => (((~g) ^ IntPart(~m) * (~g * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && IntegerQ(Simplify(~n / (~m + 1)))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (-((~b * ~e - ~a * ~f)) * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q) / (~a * ~b * ~g * ~n * (~p + 1)) + (1 / (~a * ~b * ~n * (~p + 1))) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q - 1) * Simp(~c * (~b * ~e * ~n * (~p + 1) + (~b * ~e - ~a * ~f) * (~m + 1)) + ~d * (~b * ~e * ~n * (~p + 1) + (~b * ~e - ~a * ~f) * (~m + ~n * ~q + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n], ~x) && (LtQ(~p, -1) && (GtQ(~q, 0) && Not(EqQ(~q, 1) && SimplerQ(~b * ~c - ~a * ~d, ~b * ~e - ~a * ~f))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~q * (~e + ~(f') * (~x) ^ ~n), ~x) => (-((~b * ~e - ~a * ~f)) * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~q + 1)) / (~a * ~g * ~n * (~b * ~c - ~a * ~d) * (~p + 1)) + (1 / (~a * ~n * (~b * ~c - ~a * ~d) * (~p + 1))) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q * Simp(~c * (~b * ~e - ~a * ~f) * (~m + 1) + ~e * ~n * (~b * ~c - ~a * ~d) * (~p + 1) + ~d * (~b * ~e - ~a * ~f) * (~m + ~n * (~p + ~q + 2) + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~q], ~x) && LtQ(~p, -1)
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n), ~x) => (~f * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ ~q) / (~b * ~g * (~m + ~n * (~p + ~q + 1) + 1)) + (1 / (~b * (~m + ~n * (~p + ~q + 1) + 1))) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ (~q - 1) * Simp(~c * ((~b * ~e - ~a * ~f) * (~m + 1) + ~b * ~e * ~n * (~p + ~q + 1)) + (~d * (~b * ~e - ~a * ~f) * (~m + 1) + ~f * ~n * ~q * (~b * ~c - ~a * ~d) + ~b * ~e * ~d * ~n * (~p + ~q + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (GtQ(~q, 0) && Not(EqQ(~q, 1) && SimplerQ(~e + ~f * (~x) ^ ~n, ~c + ~d * (~x) ^ ~n)))
    @apply_utils Antiderivative(((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~e + ~(f') * (~x) ^ ~n)) / (~c + ~(d') * (~x) ^ ~n), ~x) => Antiderivative(ExpandIntegrand(((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~e + ~f * (~x) ^ ~n)) / (~c + ~d * (~x) ^ ~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~q * (~e + ~(f') * (~x) ^ ~n), ~x) => ~e * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) + ((~f * (~g * ~x) ^ ~m) / (~x) ^ ~m) * Antiderivative((~x) ^ (~m + ~n) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p') * (~c + ~(d') * (~x) ^ ~(mn')) ^ ~(q') * (~e + ~(f') * (~x) ^ ~(n')) ^ ~(r'), ~x) => Antiderivative((~x) ^ (~m - ~n * ~q) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~d + ~c * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~r], ~x) && (EqQ(~mn, -(~n)) && IntegerQ(~q))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * (~c + ~(d') * (~x) ^ ~(mn')) ^ ~(q') * (~e + ~(f') * (~x) ^ ~(n')) ^ ~(r'), ~x) => Antiderivative((~x) ^ (~m + ~n * (~p + ~r)) * (~b + ~a * (~x) ^ -(~n)) ^ ~p * (~c + ~d * (~x) ^ -(~n)) ^ ~q * (~f + ~e * (~x) ^ -(~n)) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~q], ~x) && (EqQ(~mn, -(~n)) && (IntegerQ(~p) && IntegerQ(~r)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * (~c + ~(d') * (~x) ^ ~(mn')) ^ ~q * (~e + ~(f') * (~x) ^ ~(n')) ^ ~(r'), ~x) => (((~x) ^ (~n * FracPart(~q)) * (~c + ~d * (~x) ^ -(~n)) ^ FracPart(~q)) / (~d + ~c * (~x) ^ ~n) ^ FracPart(~q)) * Antiderivative((~x) ^ (~m - ~n * ~q) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~d + ~c * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~mn, -(~n)) && Not(IntegerQ(~q)))
    @apply_utils Antiderivative((~g * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p') * (~c + ~(d') * (~x) ^ ~(mn')) ^ ~(q') * (~e + ~(f') * (~x) ^ ~(n')) ^ ~(r'), ~x) => (((~g) ^ IntPart(~m) * (~g * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ -(~n)) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && EqQ(~mn, -(~n))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e + ~(f') * (~x) ^ ~n) ^ ~(r'), ~x) => Unintegrable((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p, ~q, ~r], ~x)
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * (~v) ^ ~n) ^ ~(p') * (~(c') + ~(d') * (~v) ^ ~n) ^ ~(q') * (~e + ~(f') * (~v) ^ ~n) ^ ~(r'), ~x) => ((~u) ^ ~m / (Coefficient(~v, ~x, 1) * (~v) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e + ~f * (~x) ^ ~n) ^ ~r, ~x), ~x, ~v) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~q, ~r], ~x) && LinearPairQ(~u, ~v, ~x)
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e1 + ~(f1') * (~x) ^ ~(n2')) ^ ~(r') * (~e2 + ~(f2') * (~x) ^ ~(n2')) ^ ~(r'), ~x) => Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e1 * ~e2 + ~f1 * ~f2 * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e1, ~f1, ~e2, ~f2, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~n2, ~n / 2) && (EqQ(~e2 * ~f1 + ~e1 * ~f2, 0) && (IntegerQ(~r) || GtQ(~e1, 0) && GtQ(~e2, 0))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * (~c + ~(d') * (~x) ^ ~n) ^ ~(q') * (~e1 + ~(f1') * (~x) ^ ~(n2')) ^ ~(r') * (~e2 + ~(f2') * (~x) ^ ~(n2')) ^ ~(r'), ~x) => (((~e1 + ~f1 * (~x) ^ (~n / 2)) ^ FracPart(~r) * (~e2 + ~f2 * (~x) ^ (~n / 2)) ^ FracPart(~r)) / (~e1 * ~e2 + ~f1 * ~f2 * (~x) ^ ~n) ^ FracPart(~r)) * Antiderivative((~g * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q * (~e1 * ~e2 + ~f1 * ~f2 * (~x) ^ ~n) ^ ~r, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e1, ~f1, ~e2, ~f2, ~g, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~n2, ~n / 2) && EqQ(~e2 * ~f1 + ~e1 * ~f2, 0))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.6 (g x)^m (a+b x^n)^p (c+d x^n)^q (e+f x^n)^r.jl


text:
card(stack, idx) = stack[idx]

function replace_card!(stack, r)
    stack[r.first] = r.second
    stack
end

insert_card_at_top!(stack, new_card) = push!(stack, new_card)

function remove_card_from_top!(stack)
    pop!(stack)
    stack
end

insert_card_at_bottom!(stack, new_card) = pushfirst!(stack, new_card)

function remove_card_from_bottom!(stack)
    popfirst!(stack)
    stack
end

remove_card!(stack, idx) = deleteat!(stack, idx)

check_stack_size(stack, n) = length(stack) == n

max_stars_repo_name: exercism-bot/julia

max_issues_repo_name: exercism-bot/julia

max_forks_repo_name: exercism-bot/julia


text:
# basic test of multiplying big and small together should give small

using ApproxManifoldProducts
using Test
using Manifolds
using TensorCast

##

@testset "multiply big and small" begin

##

M = SpecialEuclidean(2)
N = 100
u0 = ProductRepr([0;0.0],[1 0; 0 1.0])
ϵ = identity_element(M, u0)

X1 = [exp(M, ϵ, hat(M, ϵ, randn(3))) for i in 1:N]
X2 = [exp(M, ϵ, hat(M, ϵ, 0.1.*randn(3))) for i in 1:N]

# test get_coordinates
testval = vee(M, ϵ, log(M, ϵ, X1[1]))
@test length(testval) === 3
@test all(abs.(testval) .< 10.0)

##

p = manikde!(M, X1)
q = manikde!(M, X2)

# check new MKD have right type info cached
@test (p._u0 |> typeof) == typeof(u0)
@test (p._u0 |> typeof) == typeof(u0)

pq = manifoldProduct([p;q], M)

# check new product also has right point type info cached
@test (pq._u0 |> typeof) == typeof(u0)

##

X12_ = getPoints(pq)

X12 = AMP._pointsToMatrixCoords(pq.manifold, X12_)

@test 0.7*N < sum(abs.(X12[1,:]) .< 0.3)
@test 0.7*N < sum(abs.(X12[2,:]) .< 0.3)
@test 0.7*N < sum(abs.(X12[3,:]) .< 0.3)

##

end


text:
@testset "draw_gif_result func" begin
    @testset "using KT" begin
        K_log = 0.155  # [1/s]
        T_log = 80.5  # [s]
        u0 = 10 * (1852.0 / 3600)  # [m/s] (knot * 1852/3600)
        duration = 50  # [s]
        sampling = 1000
        time_list = range(0.0, stop = duration, length = sampling)
        Ts = 50.0
        δ_list = 10.0 * pi / 180.0 * sin.(2.0 * pi / Ts * time_list) # [rad]
        kt_results = kt_simulate(K_log, T_log, time_list, δ_list)
        r, δ = kt_results
        u = u0 * ones(Float64, length(time_list))
        v = zeros(Float64, length(time_list))
        x, y, ψ = calc_position(time_list, u, v, r)

        test_result_file_name = "test_kt.gif"
        shape = [20, 5]
        draw_gif_result(time_list, x, y, ψ, shape, test_result_file_name, fps = 5000)
        rm(test_result_file_name)
    end

end

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text:
# Testsets for Loss Functions
@testset "Cross Entropy" begin
    x = rand(Float32, 10)
    y = rand(Float32, 10)

    z = Flux.crossentropy(x, y)

    backend = nGraph.Backend()
    f = nGraph.compile(backend, Flux.crossentropy, x, y)

    @test isapprox(z, parent(f())[])
end

#####
##### Test routines for Embedding Tables
#####

unif(b) = (dims...) -> _unif(b, dims...)
function _unif(b::T, dims...) where {T}
    sampler = Uniform{T}(-b, b)
    return convert.(T, rand(sampler, dims...))
end

# function test_embedding()
#     context_size = 2
#     embedding_dim = 10
#     test_sentence = split("""
#         When forty winters shall besiege thy brow,
#         And dig deep trenches in thy beauty's field,
#         Thy youth's proud livery so gazed on now,
#         Will be a totter'd weed of small worth held:
#         Then being asked, where all thy beauty lies,
#         Where all the treasure of thy lusty days;
#         To say, within thine own deep sunken eyes,
#         Were an all-eating shame, and thriftless praise.
#         How much more praise deserv'd thy beauty's use,
#         If thou couldst answer 'This fair child of mine
#         Shall sum my count, and make my old excuse,'
#         Proving his beauty by succession thine!
#         This were to be new made when thou art old,
#         And see thy blood warm when thou feel'st it cold.
#         """)
#
#     # Build vocab and trigrams.
#     trigrams = [((test_sentence[i], test_sentence[i+1]), test_sentence[i+2]) for i in 1:length(test_sentence)-2]
#
#     vocab = Set(test_sentence)
#     word_to_idx = Dict(word => convert(Int32, i-1) for (i, word) in enumerate(vocab))
#
#     embedding_param = Flux.param(randn(Float32, embedding_dim, length(vocab)))
#     insize = context_size * embedding_dim
#
#     d1 = Dense(
#         insize,
#         128,
#         relu;
#         initW = unif(Float32(1/sqrt(insize))),
#         initb = unif(Float32(1/sqrt(insize))),
#     )
#     d2 = Dense(
#         128,
#         length(vocab);
#         initW = unif(Float32(1/sqrt(128))),
#         initb = unif(Float32(1/sqrt(128))),
#     )
#
#     f = function(x)
#         a = reshape(embedding(x, embedding_param), :)
#         # Mark the embedding parameter as an inplace update node.
#         nGraph.__inplace(embedding_param)
#         b = d1(a)
#         c = d2(b)
#         return softmax(c)
#     end
#
#     loss(a, b) = Flux.crossentropy(f(a), nGraph.onehot(b, length(vocab), 1))
#
#     # Construct inputs and output placeholders
#     inputs = zeros(Int32, context_size)
#     expected = nGraph.parameter(Int32(1))
#
#     backend = nGraph.Backend("CPU")
#     F = nGraph.compile(backend, loss, inputs, expected; optimizer = nGraph.SGD(Float32(0.01)))
#
#     # Get references to the input tensors in F
#     #inputs_view = F.inputs[1]
#     expected_view = F.inputs[2].base
#
#     context_idx = Array{Int32}(undef, context_size)
#
#     # Create a training loop
#     for epoch in 1:80
#         total_loss = 0
#         for (context, target) in shuffle(trigrams)
#             inputs .= getindex.(Ref(word_to_idx), context)
#             expected_view[] = word_to_idx[target]
#
#             this_loss = F().base
#             total_loss += this_loss[]
#         end
#         println("Epoch $epoch - Loss: $total_loss")
#     end
#
#     G = nGraph.compile(backend, f, inputs)
#
#     # Copy over the parameters from the trained model.
#     for (a, b) in zip(nGraph.getinputs(G.optimizer), nGraph.getinputs(F.optimizer))
#         @assert size(a) == size(b)
#         a .= b
#     end
#
#     return G, trigrams, word_to_idx
# end
#
# @testset "Testing Embedding Training" begin
#     fex, trigrams, word_to_idx = test_embedding()
#
#     input_view = fex.inputs[1].base
#     ncorrect = 0
#     total = length(word_to_idx)
#     for (context, target) in trigrams
#         input_view .= [word_to_idx[i] for i in context]
#         y = last(findmax(fex().base)) - 1
#         z = word_to_idx[target]
#         if y == z
#             ncorrect += 1
#         end
#     end
#     @test ncorrect / total > 0.95
# end


text:
"""
    tdc(T::SymTridiagonal)

Tridiagonal divide-and-conquer using arrowhead matrices.

COMPUTES: eigenvectors U and eigenvalues Λ of a SymTridiagonal matrix T, T = U*Diagonal(Λ)*U'.

RETURNS: Eigen(Λ, U)
"""
function tdc(T::SymTridiagonal{S}) where S

    n=length(T.dv)
    # set the tolerances for eigen
    E=Eigen(Vector{Float64}(undef,1),Array{Float64}(undef,1,1))
    τ=[1e2,1e2,1e2,1e2,1e2]
    if n==1
        E=Eigen([T.dv[1]],ones(1,1))
    elseif n==2
        E,info=eigen(SymArrow([T.dv[1]],[T.ev[1]],T.dv[2],2),τ)
    else
        k=div(n,2)
        T₁=SymTridiagonal(T.dv[1:k],T.ev[1:k-1])
        T₂=SymTridiagonal(T.dv[k+2:end],T.ev[k+2:end])
        E₁=tdc(T₁)
        E₂=tdc(T₂)
        # form arrowhead
        D=[E₁.values;E₂.values]
        z=[E₁.vectors'[:,end]*T.ev[k]; E₂.vectors'[:,1]*T.ev[k+1]]
        a=T.dv[k+1]
        A=SymArrow(D,z,a,k+1)
        E,info=eigen(A,τ)
        E.vectors[1:k,:].=E₁.vectors*E.vectors[1:k,:]
        E.vectors[k+2:end,:].=E₂.vectors*E.vectors[k+2:end,:]
    end
    return E
end

# roots

"""
    rootsah(pol::Polynomial, D::Vector,τ)

COMPUTES: the roots of polynomials with all distinct real roots.
* τ is the vector of tolerances for `eigen()`, defaults to [1e2,1e2,1e2,1e2,1e2].
* The computation is forward stable. The program uses `SymArrow` (arrowhead) companion
matrix and the corresponding `eigen()` function
* D are barycentric coordinates - elements of D must interpolate the roots of P,
for example `D=roots(polyder(pol))`

RETURNS: roots, Qout
* Qout[k] = 1 / 0 - Double was / was not used when computing k-th root
"""
function rootsah(pol::Union{Polynomial{Float32},Polynomial{Float64},Polynomial{Int32},
    Polynomial{Int64}}, D::Vector{Float64},τ::Vector{Float64}=[1e2,1e2,1e2,1e2,1e2])

    # Type is Float64
    T = Float64
    Dm=map(T,D)
    Dm=sort(Dm,rev=true)
    D=sort(D,rev=true)
    p=map(T,[pol[i] for i=0:1:length(pol)-1])
    p=p[end:-1:1]

    # Compute z of the arrowhead
    # First we compute values s=p(D) we use Horner sheme with Double
    n=length(p)-1
    pD=map(Double,p)
    DD=map(Double,Dm)
    oneD=Double(one(T))
    s=Array{Double{T}}(undef,n-1)
    for i=1:n-1
        s[i]=oneD
    end
    s=pD[1]*s
    for i=2:n+1
        r=s.*DD
        s=r.+pD[i]
    end

    # Compute t's
    t=Array{Double{T}}(undef,n-1)
    for j=1:n-1
        h=oneD
        for i=1:j-1
            g=DD[j]-DD[i]
            h=h*g
        end
        for i=j+1:n-1
            g=DD[j]-DD[i]
            h=h*g
        end
        t[j]=h
    end
    # Compute αD
    αD=pD[2]/pD[1]
    for i=1:n-1
        αD+=DD[i]
    end
    # αD=-αD
    αD=Double(-αD.hi,-αD.lo)
    #  Compute z
    zD=Array{Double{T}}(undef,n-1)
    for i=1:n-1
        #         zD[i]=sqrt((-s[i])/(t[i]*pD[1]))
        zD[i]=sqrt(Double(-s[i].hi,-s[i].lo)/(t[i]*pD[1]))
    end

    z=Array{Float64}(undef,n-1)
    for i=1:n-1
        z[i]=Float64(zD[i].lo+zD[i].hi)
    end
    α=Float64(αD.hi+αD.lo)
    A=SymArrow(D,z,α,n)
    E=zeros(Float64,n)
    Qout=zeros(Int64,n)
    for k=1:n
        E[k], Qout[k] = eigen( A,zD,αD,k,τ )
    end
    E, Qout
end

#------------------------
#--------  Functions

#-------- Inverses

"""
    inv(A::SymArrow, zD::Vector{Double},αD::Double, i, τ)

Version of `inv(A::SymArrow, i, τ)` for use with `rootsah(pol::Polynomial, D::Vector,τ)`.
"""
function inv(A::SymArrow{T},zD::Vector{Double{T}}, αD::Double{T},i::Int64,τ::Vector{Float64}=[1e3, 10*length(A.D)]) where T
    # COMPUTES: inverse of a SymArrow matrix A, inv(A-A.D[i]*I) which is again SymArrow
    # uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # τ=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used

    n=length(A.D)
    D=Array{T}(undef,n)
    z=Array{T}(undef,n)
    wz=one(T)/A.z[i]
    σ=A.D[i]
    for k=1:i-1
        D[k]=one(T)/(A.D[k]-σ)
        z[k]=-A.z[k]*D[k]*wz
    end
    for k=i+1:n
        D[k-1]=one(T)/(A.D[k]-σ)
        z[k-1]=-A.z[k]*D[k-1]*wz
    end
    a=A.a-A.D[i]

    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    nn=n-1
    for k=1:i-1
        D[k]>0.0 ? P+=A.z[k]^2*D[k] : Q+=A.z[k]^2*D[k]
    end
    for k=i:nn
        D[k]>0.0 ? P+=A.z[k+1]^2*D[k] : Q+=A.z[k+1]^2*D[k]
    end
    a<0 ? P-=a : Q-=a
    Kb=(P-Q)/abs(P+Q)
    Kz=maximum(abs,A.z)*abs(wz)
    if Kb<τ[1] ||  Kz<τ[2]
        b=(P+Q)*wz*wz
    elseif Kb==Inf
       b=0.0
    else  # recompute in Double
        Qout=1
        AD=map(Double{T},A.D)
        σd=AD[i]
        Dd=[Double{T}[AD[k]-σd for k=1:i-1];
            Double{T}[AD[k]-σd for
                            k=i+1:length(A.D)]]
        wzd=zD[i] # Double(A.z[i])
        ad=αD-σd
        Pd,Qd=map(Double{T},(0.0,0.0))
        for k=1:i-1
            Dd[k].hi>0.0 ? Pd+=zD[k]^2/Dd[k] :
            Qd+=zD[k]^2/Dd[k]
        end
        for k=i:nn
            Dd[k].hi>0.0 ? Pd+=zD[k+1]^2/Dd[k] :
            Qd+=zD[k+1]^2/Dd[k]
        end
        ad.hi<0 ?   Pd-=ad : Qd-=ad

        bd=(Pd+Qd)/(wzd*wzd)
        b=Float64(bd.hi+bd.lo)
    end

    if i<A.i
        # SymArrow([D[1:A.i-2],zero(T),D[A.i-1:end]],[z[1:A.i-2],wz,z[A.i-1:end]],b,i),Kb,Kz,Qout
        # This reduces memory allocation (4 vectors) and is 30-40% faster
        for k=n-1:-1:A.i-1
            D[k+1]=D[k]
            z[k+1]=z[k]
        end
        D[A.i-1]=zero(T)
        z[A.i-1]=wz
        SymArrow(D,z,b,i),Kb,Kz,Qout
    else
        # SymArrow([D[1:A.i-1],zero(T),D[A.i:end]],[z[1:A.i-1],wz,z[A.i:end]],b,i+1),Kb,Kz,Qout
        for k=n-1:-1:A.i
            D[k+1]=D[k]
            z[k+1]=z[k]
        end
        D[A.i]=zero(T)
        z[A.i]=wz
        SymArrow(D,z,b,i+1),Kb,Kz,Qout
    end
end # inv

"""
    inv(A::SymArrow, zD::Vector{Double},αD::Double, σ::Float64, τ)

Version of `inv(A::SymArrow, σ::Float64, τ)` for use with `rootsah(pol::Polynomial, D::Vector,τ)`.
"""
function inv(A::SymArrow{T}, zD::Vector{Double{T}}, αD::Double{T}, σ::Float64, τ::Float64=1.0e3) where T
    # COMPUTES: inverse of the shifted SymArrow A, inv(A-σ*I) which is SymDPR1
    # uses DoubleDouble to compute ρ accurately, if needed.
    # τ is tolerance, usually 1e3,  0.0 forces Double, 1e50 would never use it
    # RETURNS: SymDPR1(D,u,ρ), Kρ, Qout
    # Kρ - condition Kρ, Qout = 1 / 0 - Double was / was not used

    n=length(A.D)
    D=Array{T}(undef,n+1)
    u=Array{T}(undef,n+1)
    for k=1:A.i-1
        D[k]=one(T)/(A.D[k]-σ)
        u[k]=A.z[k]*D[k]
    end
    D[A.i]=zero(T)
    u[A.i]=-one(T)
    for k=A.i:n
        D[k+1]=one(T)/(A.D[k]-σ)
        u[k+1]=A.z[k]*D[k+1]
    end

    # compute ρ and Kρ
    #--- compute the sum in a plain loop
    a=A.a-σ
    P=zero(T)
    Q=zero(T)
    Qout=0
    for k=1:A.i-1
        D[k]>0.0 ? P+=A.z[k]^2*D[k] : Q+=A.z[k]^2*D[k]
    end
    for k=A.i:n
        D[k+1]>0.0 ? P+=A.z[k]^2*D[k+1] : Q+=!=()A.z[k]^2*D[k+1]
    end
    a<0 ? P-=a : Q-=a

    # Condition of ρ
    Kρ=(P-Q)/abs(P+Q)
    if Kρ<τ
        ρ=-1.0/(P+Q)
    else  # recompute in Double
        Qout=1
        Pd,Qd=map(Double{T},(0.0,0.0))
        AD=map(Double{T},A.D)
        σd=map(Double{T},σ)
        ad=αD-σd
        for k=1:A.i-1
            D[k]>0.0 ? Pd+=zD[k]^2/(AD[k]-σd) : Qd+=zD[k]^2/(AD[k]-σd)
        end
        for k=A.i:n
            D[k+1]>0.0 ? Pd+=zD[k]^2/(AD[k]-σd) : Qd+=zD[k]^2/(A.D[k]-σd)
        end
        ad.hi+ad.lo<0 ? Pd=Pd-ad : Qd=Qd-ad
        r=map(Double{T},1.0)/(Pd+Qd)
        ρ=Float64(-(r.hi+r.lo))
    end

    # returns the following
    SymDPR1(D,u,ρ), Kρ, Qout
end # inv

"""
    inv(A::SymArrow, zD::Vector{Double},αD::Double, σ::Double)

Version of `inv(A::SymArrow, σ::Double)` for use with `rootsah(pol::Polynomial, D::Vector,τ)`.
"""
function inv(A::SymArrow{T}, zD::Vector{Double{T}}, αD::Double{T}, σ::Double{Float64}) where T
    # COMPUTES: inverse of the σed SymArrow A, inv(A-σ*I), which is a SymDPR1.
    # Here σ is Double so it uses Double to compute everything.
    # RETURNS: SymDPR1(D1,u1,ρ), Qout
    # Qout = 1 on exit meaning Double was used

    n=length(A.D)
    D=Array{Double{T}}(undef,n+1)
    u=Array{Double{T}}(undef,n+1)
    σd=map(Double{T},σ)
    for k=1:n
        D[k]=map(Double{T},A.D[k])-σd
    end
    u[1:n]=zD
    a=αD-σd
    i=A.i

    oned=Double(one(T))
    zerod=Double(zero(T))
    for k=1:i-1
        u[k]=u[k]/D[k]
        D[k]=oned/D[k]
    end
    for k=i:n
        u[k+1]=u[k]/D[k]
        D[k+1]=oned/D[k+1]
    end
    D[i]=zerod
    u[i]=-oned

    # compute ρ and Kρ
    #--- compute the sum in a plain loop
    P,Q=zerod,zerod
    Qout=1
    for k=1:i-1
        D[k].hi > 0.0 ? P+=zD[k]*u[k] : Q+=zD[k]*u[k]
    end
    for k=i:n
        D[k+1].hi >0.0 ? P+=zD[k]*u[k+1] : Q+=zD[k]*u[k+1]
    end
    a.hi+a.lo<0.0  ? P-=a : Q-=a
    r=oned/(P+Q)
    ρ=Float64(-(r.hi+r.lo))

    # returns the following
    D₁=Array{T}(undef,n+1)
    u₁=Array{T}(undef,n+1)
    for k=1:n+1
        D₁[k]=Float64(D[k].hi+D[k].lo)
    end
    for k=1:n+1
        u₁[k]=Float64(u[k].hi+u[k].lo)
    end
    SymDPR1(D₁,u₁,ρ), Qout
end # inv

"""
    eigen(A::SymArrow, zD::Vector{Double}, αD::Double, k, τ)

Version of `eigen(A::SymArrow, k, τ)` for use with `rootsah(pol::Polynomial, D::Vector,τ)`.
"""
function  eigen( A::SymArrow{T},zD::Vector{Double{T}},αD::Double{T},k::Int64,
    τ::Vector{Float64}=[1e3,10.0*length(A.D),1e3,1e3,1e3]) where T
    # COMPUTES: k-th eigenvalue of an ordered irreducible SymArrow
    # A = [Diagonal(D) z; z' α]
    # Specially designed to be used in the polynomial rootfinder rootsah !!!!
    # τ=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: λ
    # λ - k-th eigenvalue
    # Kb, Kz, Kν, Kρ - condition numbers
    # Qout = 1 / 0 - Double was / was not used
    # If Qout>0, quad precision was used
    # i was the shift index

    # Set the dimension
    n = length(A.D) + 1
    # Set all conditions initially to zero
    Kb,Kz,Kν,Kρ=0.0,0.0,0.0,0.0
    Qout=0
    v=zeros(T,n)
    # Kz is former kappa_nu

    # Determine the shift σ, the shift index i, and whether λ
    # is on the left or the right side of the nearest pole
    # Exterior eigenvalues (k = 1 or k = n):
    if k == 1
        σ,i,side = A.D[1],1,'R'
    elseif k==n
        σ,i,side = A.D[n-1],n-1,'L'
    else
        # Interior eigenvalues (k in (2,...n-1) ):
        # We need to compute this in Double
        Dd=map(Double,A.D)
        middle=(Dd[k-1]-Dd[k])/Double(2.0)
        Fmiddle=Double(zero(T))
        for l=1:n-1
            Fmiddle+=zD[l]^2/((Dd[l]-Dd[k])-middle)
        end
        Fmiddle=((αD-Dd[k])-middle)-Fmiddle
        σ,i,side = Fmiddle.hi+Fmiddle.lo < 0.0 ? (A.D[k],k,'R') : (A.D[k-1],k-1,'L')
    end

    # Compute the inverse of the shifted matrix, A_i^(-1), Kb and Kz
    Ainv,Kb,Kz,Qout = inv(A,zD,αD,i)

    # Compute the eigenvalue of the inverse shifted matrix
    ν = bisect( Ainv,side )

    #  ν=fastzero([invD1; 0; invD2], [w1;wz;w2], b, side); # implement later
    #  [ν-ν₁]/ν, [ν-nueva]/ν, pause  # and compare
    if abs(ν)==Inf
        # this is nonstandard
        # Deflation in aheig (ν=Inf)
        v[i]=1.0
        λ=σ
    else
        # standard case, full computation
        # ν₁ is the F- or 1-norm of the inverse of the shifted matrix
        # ν₁0=maximum([sum(abs,Ainv.z)+abs(Ainv.a), maximum(abs.(Ainv.D)+abs.(Ainv.z))])
        ν₁=0.0
        for l=1:n-1
            ν₁=max(ν₁,abs(Ainv.D[l])+abs(Ainv.z[l]))
        end
        ν₁=max(sum(abs,Ainv.z)+abs(Ainv.a), ν₁)
        Kν=ν₁/abs(ν)
        if Kν<τ[3]
            # Accuracy is fine, compute the eigenvector
            μ = 1.0/ν
            # v=[ A.z./((A.D-σ)-μ);-1.0]
            for l=1:n-1
            	v[l] = A.z[l]/((A.D[l]-σ)-μ)
            end
            v[n]=-1.0
            λ=μ+σ
            normalize!(v)
        else
            # Remedies according to Remark 3 - we shift between original
            # eigenvalues and compute DPR1 matrix
            # 1/ν₁+σ, 1/ν+σ
            # println(k, " Remedy 3 ")
            ν = side=='R' ? abs(ν) : -abs(ν)
            ν₁=-sign(ν)*ν₁
            σ₁=(ν₁+ν)/(2.0*ν*ν₁)+σ
            if count(isequal(σ₁),A.D)>0 # we came back with a pole
                # recompute σᵥ more accurately according with dekker
                σᵥ=(Double(ν₁)+Double(ν))/(Double(2.0)*Double(ν)*Double(ν₁))+Double(σ)
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv,Qout₁=inv(A,σᵥ) # Ainv is Float64
                ν₁=bisect(Ainv,side)
                μ₁ = 1.0/ν₁
                D0=map(Double,A.D)-σᵥ
                D1=[D0[l].hi+D0[l].lo for l=1:length(D0)]
                if count(isequal(μ₁),D1)>0
                    ind=findall(isequal(μ₁),D1)
                    v=zeros(n)
                    v[ind]=1.0
                else
                    v=[ A.z./(D1-μ₁);-1.0]
                end
                # Shift the eigenvalue back in Double
                lam = Double(1.0)/Double(ν₁)+σᵥ
                # Return this
                λ = lam.hi+lam.lo
                normalize!(v)
            else
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv, Kρ,Qout₁=inv(A,σ₁,τ[4]) # Ainv is Float64
                # Compute the eigenvalue by bisect for DPR1
	            # Note: instead of bisection could use dlaed4 (need a wrapper) but
	            # it is not faster. There norm(u)==1
                ν₁= bisect(Ainv,side)
                μ₁=1.0/ν₁
                # standard v
                v=[ A.z./((A.D-σ₁)-μ₁);-1.0]
                # Return this - shift the eigenvalue back and normalize the vector
                λ = μ₁+σ₁
                normalize!(v)
            end
            Qout==1 && (Qout=Qout+2)
        end

        # Remedy according to Remark 1 - we recompute the the eigenvalue
        # near zero from the inverse of the original matrix (a DPR1 matrix).
        if (abs(A.D[i])+abs(1.0/ν))/abs(λ)>τ[5]
            if k==1 && A.D[1]<0.0 || k==n && A.D[n-1]>0.0 || i<n-1 && side=='L' && sign(A.D[i])+sign(A.D[i+1])==0 ||
                side=='R' && sign(A.D[i])+sign(A.D[i-1])==0
                # println(k, " Remedy 1 ")
                # Compute the inverse of the original arrowhead (DPR1)
                Ainv,Kρ,Qout₁ = inv(A,0.0,τ[4]) # Ainv is Float64
                Qout==1 && (Qout=Qout+4)
                if abs(Ainv.r)==Inf
                    λ=0.0
                else
                    # Here we do not need bisection. We compute the Rayleigh
                    # quotient by using already computed vectors which is
                    # componentwise accurate
                    ν₁=sum(v.^2 .*Ainv.D)+Ainv.r*sum(v.*Ainv.u)^2;
                    λ=1.0/ν₁
                end
            end
        end
    end

    # Return this
    λ, Qout
end # eigen (k)


text:
function test_handler(session, req)
    global test_observable
    test_observable = Observable(Dict{String, Any}())
    test_session = session

    s1 = JSServe.Slider(1:100)
    s2 = JSServe.Slider(1:100)
    b = JSServe.Button("hi"; dataTestId="hi_button")

    clicks = Observable(0)
    on(b) do click
        clicks[] = clicks[] + 1
    end
    clicks_div = DOM.div(clicks, dataTestId="button_clicks")
    t = JSServe.TextField("Write!")

    linkjs(session, s1.value, s2.value)

    onjs(session, s1.value, js"""function (v){
        var updated = JSServe.update_obs($(test_observable), {onjs: v});
        console.log(updated);
    }""")

    on(t) do value
        test_observable[] = Dict{String, Any}("textfield" => value)
    end

    on(b) do value
        test_observable[] = Dict{String, Any}("button" => value)
    end

    textresult = DOM.div(t.value; dataTestId="text_result")
    sliderresult = DOM.div(s1.value; dataTestId="slider_result")

    number_input = JSServe.NumberInput(66.0)
    number_result = DOM.div(number_input.value, dataTestId="number_result")

    dom = md"""
    # IS THIS REAL?

    My first slider: $(s1)

    My second slider: $(s2)

    Test: $(sliderresult)

    The BUTTON: $(b)

    $(clicks_div)

    # Number Input

    $(number_input)

    $(number_result)

    Type something for the list: $(t)

    some list $(textresult)
    """
    return DOM.div(dom)
end

function test_current_session(app)
    dom = children(app.dom)[1]
    @test evaljs(app, js"document.querySelectorAll('input[type=\"button\"]').length") == 1
    @test evaljs(app, js"document.querySelectorAll('input[type=\"range\"]').length") == 2
    @test evaljs(app, js"document.querySelectorAll('input[type=\"button\"]').length") == 1
    @test evaljs(app, js"document.querySelectorAll('input[type=\"range\"]').length") == 2

    @testset "button" begin
        # It's in the dom!
        @test evaljs(app, js"document.querySelectorAll('input[type=\"button\"]').length") == 1
        bquery = query_testid("hi_button")
        # Spam the button press on the JS side a bit, to make sure we're not loosing events!
        for i in 1:100
            val = test_value(app, js"$(bquery).click()")
            @test val["button"] == true
        end
        button = query_testid("button_clicks")
        @test evaljs(app, js"$(button).innerText") == "100"
        button = dom.content[5].content[2]
        @test button.content[] == "hi"
        button.content[] = "new name"
        bquery = query_testid("hi_button")
        @test evaljs(app, js"$(bquery).value") == "new name"
        # button press from Julia
        val = test_value(app, ()-> (button.value[] = true))
        @test val["button"] == true
    end

    @testset "textfield" begin
        @test evaljs(app, js"document.querySelectorAll('input[type=\"textfield\"]').length") == 1
        # @test evaljs(app, js"document.querySelector('input[type=\"textfield\"]').value") == "Write!"
        # Spam the button press a bit!
        text_obs = children(dom.content[end].content[2])[1]
        textfield = dom.content[end-1].content[2]
        text_result = query_testid("text_result")
        @test evaljs(app, js"$(text_result).innerText") == "Write!"
        @testset "setting value from js" begin
            for i in 1:10
                str = randstring(10)
                do_input = js"""
                    var tfield = document.querySelector('input[type=\"textfield\"]');
                    tfield.value = $(str);
                    tfield.onchange();
                """
                val = test_value(app, do_input)
                @test val["textfield"] == str
                @test text_obs[] == str
                evaljs(app, js"$(text_result).innerText") == str
                @test textfield[] == str
            end
        end
        @testset "setting value from julia" begin
            for i in 1:10
                str = randstring(10)
                val = test_value(app, ()-> textfield[] = str)
                @test val["textfield"] == str
                @test text_obs[] == str
                evaljs(app, js"$(text_result).innerText") == str
                @test textfield[] == str
            end
        end
    end

    @testset "slider" begin
        # We test with JSCall this time, to test it as well ;)
        slider1 = dom.content[2].content[2]
        slider2 = dom.content[3].content[2]
        slider1_js = js"document.querySelectorAll('input[type=\"range\"]')[0]"
        slider2_js = js"document.querySelectorAll('input[type=\"range\"]')[1]"
        sliderresult = query_testid("slider_result")

        @testset "set via julia" begin
            for i in 1:2:100
                slider1[] = i
                @test evaljs(app, js"$(slider1_js).value") == "$i"
                # Test linkjs
                @test evaljs(app, js"$(slider2_js).value") == "$i"
                @test slider2[] == i
                evaljs(app, js"$(sliderresult).innerText") == "$i"
            end
        end
    end
end

global test_session = nothing
global dom = nothing
inline_display = JSServe.App() do session, req
    global test_session = session
    global dom = test_handler(session, req)
    return DOM.div(JSTest, dom)
end;
electron_disp = electrondisplay(inline_display);
app = TestSession(URI("http://localhost:8555/show"),
                  JSServe.GLOBAL_SERVER[], electron_disp, test_session)
app.dom = dom

@testset "electron inline display" begin
    test_current_session(app)
end
close(app)

@testset "webio mime" begin
    ENV["JULIA_WEBIO_BASEURL"] = "https://google.de/"
    JSServe.__init__()
    @test JSServe.JSSERVE_CONFIGURATION.external_url[] == "https://google.de"
    @test JSServe.JSSERVE_CONFIGURATION.content_delivery_url[] == "https://google.de"
    html_webio = sprint(io-> show(io, MIME"application/vnd.jsserve.application+html"(), inline_display))
    #@test occursin("proxy_url = 'https://google.de';", html_webio)
    # @test JSServe.url("/test") == "https://google.de/test"
    JSServe.JSSERVE_CONFIGURATION.external_url[] = ""
    JSServe.JSSERVE_CONFIGURATION.content_delivery_url[] = ""
    # @test JSServe.url("/test") == "/test" # back to relative urls
    html_webio = sprint(io-> show(io, MIME"application/vnd.jsserve.application+html"(), inline_display))
    # We open the display server with the above TestSession
    # TODO electrontests should do this!
    check_and_close_display()
end

@testset "Electron standalone" begin
    testsession(test_handler, port=8555) do app
        test_current_session(app)
    end
end

@testset "markdown" begin
    function test_handler(session, req)
        global dom, test_session, test_observable
        test_session = session
        dom = md"""
        # More MD

        [Github-flavored Markdown info page](http://github.github.com/github-flavored-markdown/)

        [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)

        Lalala
        ======

        Alt-H2
        ------

        Emphasis, aka italics, with *asterisks* or _underscores_.

        Strong emphasis, aka bold, with **asterisks** or __underscores__.

        Combined emphasis with **asterisks and _underscores_**.

        Strikethrough uses two tildes. ~~Scratch this.~~

        1. First ordered list item
        2. Another item
            * Unordered sub-list.
        1. Actual numbers don't matter, just that it's a number
            1. Ordered sub-list
        4. And another item.

        * Unordered list can use asterisks
        - Or minuses
        + Or pluses

        Inline `code` has `back-ticks around` it.

        ```julia
        test("haha")
        ```

        | Tables        | Are           | Cool  |
        | ------------- |:-------------:| -----:|
        | col 3 is      | right-aligned | $1600 |
        | col 2 is      | centered      |   $12 |
        | zebra stripes | are neat      |    $1 |

        > Blockquotes are very handy in email to emulate reply text.
        > This line is part of the same quote.

        Three or more...

        ---

        Hyphens[^1]

        ***

        Asterisks

        ___

        Underscores

        [^1]: This is the first footnote.
        """
        return DOM.div(JSServe.MarkdownCSS, dom)
    end

    testsession(test_handler, port=8555) do app
        # Lets not be too porcelainy about this ...
        md_js_dom = js"document.getElementById('application-dom').children[0]"
        @test evaljs(app, js"$(md_js_dom).children.length") == 1
        md_children = js"$(md_js_dom).children[0].children[0].children"
        @test evaljs(app, js"$(md_children).length") == 23
        @test occursin("This is the first footnote.", evaljs(app, js"$(md_children)[22].innerText"))
        @test evaljs(app, js"$(md_children)[2].children[0].children[0].tagName") == "IMG"
    end
end

@testset "range slider" begin
    function test_handler(session, req)
        rslider = JSServe.RangeSlider(1:100; value=[10, 80])
        start = map(first, rslider)
        stop = map(last, rslider)
        return DOM.div(rslider, start, stop, id="rslider")
    end
    testsession(test_handler, port=8555) do app
        # Lets not be too porcelainy about this ...
        rslider = getfield(app.dom, :children)[1]
        @test rslider[] == [10, 80]
        rslider_html = js"document.getElementById('rslider')"
        @test evaljs(app, js"$(rslider_html).children.length") == 3
        @test evaljs(app, js"$(rslider_html).children[1].innerText") == "10"
        @test evaljs(app, js"$(rslider_html).children[2].innerText") == "80"
        rslider[] = [20, 70]
    end
end


text:
@testset "sf -> node" begin
    function normal_model(t, x)
        loc = sample(t, "loc", Normal())
        log_scale = sample(t, "log_scale", Normal())
        observe(t, "data", Normal(loc, exp(log_scale)), x)
    end

    sf = TypedSF(normal_model, Maybe{Float64}, Float64)
    input = 3.0
    csf = compile(sf, input, ["log_scale", "loc"])

    t = trace()
    sample(t, "test", csf)
    @info t
    @test abs(t["test"].value["loc"] - input) < 0.1
end


text:
module Bloom
	import Base.contains, Base.show
	export BloomFilter, add!
	include("bloom-filter.jl")
end


text:
using RunTests
using Base.Test

@testmodule CommonsTest begin
  using MLLibrary
  function test_backsolve()

    A= MLLibrary.DataManipulation.DataTransform.convertToFloat([1 2 3; 0 4 6;0 0 9])
    B= MLLibrary.DataManipulation.DataTransform.convertToFloat([10;24;18])

    answer = [-2. 3. 2.]'
    @test MLLibrary.LinearAlgebra.Commons.backsolve(A,B)==answer

    A= MLLibrary.DataManipulation.DataTransform.convertToFloat([1 2 3; 0 4 6;0 0 9;0 0 0;0 0 0])
    @test MLLibrary.LinearAlgebra.Commons.backsolve(A,B)==answer

  end
end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
"""
    setup_worker(host, port[, cookie];
        threadlevel=:serialized, stdout_to_master=true, stderr_to_master=true)

This is the entrypoint for MPI workers using TCP transport.

1. it connects to the socket on master
2. sends the process rank and size
3. hands over control via [`Distributed.start_worker`](https://docs.julialang.org/en/v1/stdlib/Distributed/#Distributed.start_worker)
"""
function setup_worker(host::Union{Integer, String}, port::Integer, cookie::Union{String, Symbol, Nothing}=nothing;
        threadlevel=:serialized, stdout_to_master=true, stderr_to_master=true)
    # Connect to the manager
    ip = host isa Integer ? IPv4(host) : parse(IPAddr, host)

    io = connect(ip, port)
    wait_connected(io)
    stdout_to_master && redirect_stdout(io)
    stderr_to_master && redirect_stderr(io)

    MPI.Initialized() || MPI.Init(;threadlevel=threadlevel)
    rank = MPI.Comm_rank(MPI.COMM_WORLD)
    nprocs = MPI.Comm_size(MPI.COMM_WORLD)
    Serialization.serialize(io, rank)
    Serialization.serialize(io, nprocs)

    # Hand over control to Base
    if isnothing(cookie)
        Distributed.start_worker(io)
    else
        if isa(cookie, Symbol)
            cookie = string(cookie)[8:end] # strip the leading "cookie_"
        end
        Distributed.start_worker(io, cookie)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    filter(pred, data)

Retain all locations in geospatial `data` according to
a predicate function `pred`. A predicate function takes
table rows as input, e.g. `pred(r) = r.state == "CA"`.
"""
function Base.filter(pred, data::D) where {D<:Data}
  𝒯 = values(data)
  𝒟 = domain(data)

  # row table view
  ctor = Tables.materializer(𝒯)
  rows = Tables.rows(𝒯)

  # indices to retain
  inds = findall(pred, rows)

  # return point set
  dom = view(𝒟, inds)
  tab = ctor(rows[inds])
  dat = Dict(paramdim(dom) => tab)

  constructor(D)(dom, dat)
end


text:
include("default_params.jl")
include("model_def.jl")

include("../helper_functions.jl")

orders = orders[end]

# initial condition
d0_map_point_mass(dq) = (interior_point_mass(2.0,4,dq).coeffs)

# build the generators
models = make_approximations(orders,approx_types,d0_map_point_mass,(args...)->args[1])

# points at which to evaluate the first return cdf
x_vals = range(0.0,1.6;length=n_err_evals)

for k1 in orders
    for k2 in keys(models[k1])
        if k2!="dg2"
            @show k1
            @show k2
            println()
            # evaluate psi for all approximations
            _ffq = ffq(models[k1][k2]["dq"].mesh,Δ(models[k1][k2]["dq"],1))
            D = InOutGenerator(_ffq,0.0)
            psi = build_psi(D;err=9e-9)
            JSON.write(
                (@__DIR__)*"/data/psi/order_"*string(k1)*"_model_"*string(k2)*".json",
                JSON.json(psi),
            )

            # evaluate first-return distribution
            first_return_coeffs = zeros(size(models[k1][k2]["coeffs_0"]))
            first_return_coeffs[index(_ffq,Minus)] =
                transpose(models[k1][k2]["coeffs_0"][index(_ffq,Plus)])*psi
            first_return_dist = SFMDistribution(first_return_coeffs,_ffq.dq)
            cdf_fun = cdf(first_return_dist)
            first_return_cdf_phase_2 = cdf_fun.(x_vals,2)
            first_return_cdf_phase_4 = cdf_fun.(x_vals,4)
            CSV.write(
                (@__DIR__)*"/data/first_return_cdf_approximations/order_"*string(k1)*"_model_"*string(k2)*".csv",
                DataFrame(
                    x = x_vals,
                    phase_2 = first_return_cdf_phase_2,
                    phase_4 = first_return_cdf_phase_4,
                ),
            )
        end
    end
end


text:
using ThreeBodyDecay
using Test
using Parameters

function QTB_mismatch_factor(dc)
  k = dc.k; i,j = ij_from_k(k);
  @unpack Hij, HRk = dc
  @unpack two_s, tbs = dc
  @unpack two_js = tbs
  #
  avHHsq =
      sum((tbs.two_js[4]!=(two_τ-two_λs[k]) ? 0.0 : 1.0) *
          abs2(amplitude(Hij, two_τ, two_λs[k])) *
          abs2(amplitude(Hij, two_λs[i], two_λs[j]))
          for two_λs in itr(tbs.two_js),
              two_τ in -dc.two_s:2:dc.two_s)
  return avHHsq
end

@testset "QaudGK vs Cuba within 0.01%" begin
  ms = ThreeBodyMasses(m0=6.0,m1=1.0,m2=1.5,m3=2.0)
  #
  v_qtv = RhoQTB(ms.m0^2,σ->1.0, ms^2)
  v_fll = real(three_body_phase_space_integral(σs->1.0, ms))
  #
  @test (v_qtv-v_fll)/v_fll < 1e-5
end

@testset "QaudGK vs Cuba for P-wave within 0.01%" begin
  tbs = ThreeBodySystem(1.0,1.5,2.0, m0=6.0)
  #
  ch = DecayChainLS(1,σ->1.0; two_s=2, tbs=tbs, parity='-', Ps=['+','+','+','+'])
  sp = summed_over_polarization((σs,two_λ)->abs2(amplitude(σs,two_λ,ch)), tbs.two_js)
  #
  v_qtv = RhoQTB(tbs.ms.m0^2,σ->1.0, tbs.ms^2)
  v_qtv *= QTB_mismatch_factor(ch) # helicity couplings
  # @show v_qtv
  v_fll = real(three_body_phase_space_integral(sp, tbs.ms))
  #
  @test abs(v_qtv-v_fll)/v_fll < 1e-5
end

tbs = ThreeBodySystem(1.0,1.5,2.0; m0=6.0, two_js=ThreeBodySpins(0, 1, 0; two_h0=1))

@testset "QaudGK vs Cuba for half-integer spin within 0.01%" begin
  #
  ch = DecayChainLS(1,σ->1.0; two_s=1, tbs=tbs, parity='-', Ps=['-','+','-','+'])
  sp = summed_over_polarization((σs,two_λ)->abs2(amplitude(σs,two_λ,ch)), tbs.two_js)
  #
  v_qtv = RhoQTB(tbs.ms.m0^2,σ->1.0, tbs.ms^2)
  v_qtv *= QTB_mismatch_factor(ch) # helicity couplings
  # @show v_qtv
  v_fll = real(three_body_phase_space_integral(sp, tbs.ms))
  # @show v_fll
  #
  @test abs(v_qtv-v_fll)/v_fll < 1e-5
end

@testset "QaudGK vs Cuba for half-integer spin and P-wave within 0.01%" begin
  #
  ch = DecayChainLS(1,σ->1.0; two_s=1, tbs=tbs, parity='+', Ps=['-','+','-','+'])
  sp = summed_over_polarization((σs,two_λ)->abs2(amplitude(σs,two_λ,ch)), tbs.two_js)
  # sp = σs->abs2(amplitude(σs,[0,0,0,0],ch))
  #
  v_qtv = RhoQTB(tbs.ms.m0^2,σ->1.0, tbs.ms^2)
  v_qtv *= QTB_mismatch_factor(ch) # helicity couplings
  # @show v_qtv
  v_fll = real(three_body_phase_space_integral(sp, tbs.ms))
  # @show v_fll
  #
  @test abs(v_qtv-v_fll)/v_fll < 1e-5
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay