text:
##################################################################
##  Chapter 1: Loading the damn thing.
#
# File for parsing the RubiRules and dumping them into a pretty printed JSON
# in your Rubin.jl/src dir.
#
# Here are the commands to download the artifacts:
# RUN THEM!
# using Pkg
# Pkg.instantiate()
##################################################################

using Rubin
using Test
using Artifacts
using FileTrees
using JSON3
using StructTypes
using Logging
using SymPy
using PyCall
using CSV
using DelimitedFiles
using DataFrames

# TODO mark as constants
rubi = artifact"Rubi"
rubitests = artifact"RubiTests"
rubidir = FileTree(rubi);
abstract type AbstractRubiParser end
struct RubiRules <: AbstractRubiParser end
struct RubiTests <: AbstractRubiParser end

## Design goals:
#1. Read Rules into a huge JSON array (with FileTrees.jl)
#2. Add metadata, facts, number, filename
#3. Save to a single file

## Steps:
# 1. Read artifact dir into a FileTree
# 2. Come up with a function to parse the files
# INPUT examples:
# (* ::Subsection::Closed:: *)
# (* 1.1.1.1 (a+b x)^m *)
# Int[1/x_, x_Symbol] := Log[x]
# Int[x_^m_., x_Symbol] := x^(m + 1)/(m + 1) /; FreeQ[m, x] && NeQ[m, -1]
# (Int[...]) := (...) /; (...)
# ^ lhs          ^rhs     ^facts
#
# BUT we can also have
# (* Int[Sqrt[a_.+b_.*x_]*(A_.+B_.*x_)/(Sqrt[c_.+d_.*x_]*Sqrt[e_.+f_.*x_ ]*Sqrt[g_.+h_.*x_]),x_Symbol] :=  B*Sqrt[a+b*x]*Sqrt[e+f*x]*Sqrt[g+h*x]/(f*h*Sqrt[c+d*x]) - B*(b*g-a*h)/(2*f*h)*Int[Sqrt[e+f*x]/(Sqrt[a+b*x]*Sqrt[c+d*x]*Sqrt[g+ h*x]),x] + B*(d*e-c*f)*(d*g-c*h)/(2*d*f*h)*Int[Sqrt[a+b*x]/((c+d*x)^(3/2)*Sqrt[ e+f*x]*Sqrt[g+h*x]),x] /; FreeQ[{a,b,c,d,e,f,g,h,A,B},x] &&  EqQ[2*A*d*f-B*(d*e+c*f),0] *)
# (* Int[...] := (...) ;/ (...) *)
# ^comment syntax around lhs := rhs ;/ facts *) closing comment syntax
# Which means we have 2 grammars:
# Case 1:
# (^Int.?) := (.+) /; (.+)
# Case 2:
# (* (Int.?) := (.+)(|\/; (.+))
#
## SOLUTION:
# With the help of Regex101.com:
# (Int.+) := (.+)(?: \/; \*\))
# ^glob lhs,
#             ^glob rhs
#                  ^glob facts, if any exist, ignore if it doesn't
# Let's now write some tests for the regex.
regextest = [
"""(* ::Subsection::Closed:: *)""",
"""(* 1.1.1.1 (a+b x)^m *)""",
"""Int[1/x_, x_Symbol] := Log[x]""",
"""(* Int[Sqrt[a_.+b_.*x_]*(A_.+B_.*x_)/(Sqrt[c_.+d_.*x_]*Sqrt[e_.+f_.*x_ ]*Sqrt[g_.+h_.*x_]),x_Symbol] :=  B*Sqrt[a+b*x]*Sqrt[e+f*x]*Sqrt[g+h*x]/(f*h*Sqrt[c+d*x]) - B*(b*g-a*h)/(3*f*h)*Int[Sqrt[e+f*x]/(Sqrt[a+b*x]*Sqrt[c+d*x]*Sqrt[g+ h*x]),x] + B*(d*e-c*f)*(d*g-c*h)/(2*d*f*h)*Int[Sqrt[a+b*x]/((c+d*x)^(3/2)*Sqrt[ e+f*x]*Sqrt[g+h*x]),x] /; FreeQ[{a,b,c,d,e,f,g,h,A,B},x] &&  EqQ[2*A*d*f-B*(d*e+c*f),0] *)""",
"""(* Int[Sqrt[a_.+b_.*x_]*(A_.+B_.*x_)/(Sqrt[c_.+d_.*x_]*Sqrt[e_.+f_.*x_ ]*Sqrt[g_.+h_.*x_]),x_Symbol] :=  B*Sqrt[a+b*x]*Sqrt[e+f*x]*Sqrt[g+h*x]/(f*h*Sqrt[c+d*x]) - B*(b*g-a*h)/(2*f*h)*Int[Sqrt[e+f*x]/(Sqrt[a+b*x]*Sqrt[c+d*x]*Sqrt[g+ h*x]),x] + B*(d*e-c*f)*(d*g-c*h)/(2*d*f*h)*Int[Sqrt[a+b*x]/((c+d*x)^(3/2)*Sqrt[ e+f*x]*Sqrt[g+h*x]),x] /; FreeQ[{a,b,c,d,e,f,g,h,A,B},x] &&  EqQ[2*A*d*f-B*(d*e+c*f),0] *)""",
"""(* Int[1/(a_+b_.*x_^5),x_Symbol] := With[{r=Numerator[Rt[a/b,5]],  s=Denominator[Rt[a/b,5]]}, r/(5*a)*Int[1/(r+s*x),x] + 2*r/(5*a)*Int[(r-1/4*(1-Sqrt[5])*s*x)/(r^2-1/2*(1-Sqrt[5])*r*s*x+s^ 2*x^2),x] + 2*r/(5*a)*Int[(r-1/4*(1+Sqrt[5])*s*x)/(r^2-1/2*(1+Sqrt[5])*r*s*x+s^ 2*x^2),x]] /; FreeQ[{a,b},x] && PosQ[a/b] *)""",
"""(* Int[1/Sqrt[a_+b_.*x_^3],x_Symbol] := With[{q=Rt[b/a,3]}, -Sqrt[2]*(1+Sqrt[3])*(1+Sqrt[3]+q*x)^2*Sqrt[(1+q^3*x^3)/(1+Sqrt[3]+ q*x)^4]/(3^(1/4)*q*Sqrt[a+b*x^3])* EllipticF[ArcSin[(-1+Sqrt[3]-q*x)/(1+Sqrt[3]+q*x)],-7-4*Sqrt[3]]]  /; FreeQ[{a,b},x] && PosQ[a] *)""",
"""(* Int[Sqrt[c_+d_.*x_^2]/((a_+b_.*x_^2)*Sqrt[e_+f_.*x_^2]),x_Symbol] :=   Sqrt[c+d*x^2]*Sqrt[c*(e+f*x^2)/(e*(c+d*x^2))]/(a*Rt[d/c,2]*Sqrt[e+f* x^2])* EllipticPi[1-b*c/(a*d),ArcTan[Rt[d/c,2]*x],1-c*f/(d*e)] /;  FreeQ[{a,b,c,d,e,f},x] && PosQ[d/c] *)"""];

## Where the 🌠 magic happens 🌠
intrulesregexfacts = r"(Int.+) := (.+) \/; (.+)(?: \*\))";
intrulesregexsimple = r"^(Int.+) := (.+)";
commenttest = """(* Int[Sqrt[a_.+b_.*x_]*(A_.+B_.*x_)/(Sqrt[c_.+d_.*x_]*Sqrt[e_.+f_.*x_ ]*Sqrt[g_.+h_.*x_]),x_Symbol] :=  B*Sqrt[a+b*x]*Sqrt[e+f*x]*Sqrt[g+h*x]/(f*h*Sqrt[c+d*x]) - B*(b*g-a*h)/(2*f*h)*Int[Sqrt[e+f*x]/(Sqrt[a+b*x]*Sqrt[c+d*x]*Sqrt[g+ h*x]),x] + B*(d*e-c*f)*(d*g-c*h)/(2*d*f*h)*Int[Sqrt[a+b*x]/((c+d*x)^(3/2)*Sqrt[ e+f*x]*Sqrt[g+h*x]),x] /; FreeQ[{a,b,c,d,e,f,g,h,A,B},x] &&  EqQ[2*A*d*f-B*(d*e+c*f),0] *)""";

## Regex tests
# Gotta make sure we are slurping up the things we care about, without too much cruft.
vregex = match.(intrulesregexsimple, regextest);
rulesregex = match.(intrulesregexfacts, regextest);
@test isnothing(vregex[1])
@test isnothing(vregex[2])
@test vregex[3].captures[1] == "Int[1/x_, x_Symbol]"
@test vregex[3].captures[2] == "Log[x]"
@test rulesregex[end].captures[3] == " FreeQ[{a,b,c,d,e,f},x] && PosQ[d/c]"

"""
# Structs for parsing
- Useful in conjunction with StructTypes.jl and dumping everything into a JSON
- Fields:
	pathname: String
	filename: String
	rulenumber: Int
	commented: Bool
	lhs: String
	rhs: String
	givens: String (not atomize)
"""
Base.@kwdef struct IntRuleCapture
	pathname::String = ""
	filename::String = ""
	rulenumber::Int = ""
	comment::Bool = true
	rhs::String = ""
	lhs::String = ""
	givens::String	 = ""
end

## Auxiliary functions for later
lhs(cap::RegexMatch) = cap.captures[1]
rhs(cap::RegexMatch) = cap.captures[2]
givens(cap::RegexMatch) = cap.captures[3]
iscommented(cap::RegexMatch) = startswith(cap.captures[1], '(')

# Input: A file
# Output: An array of structs like IntRuleCapture
# TODO: Figure out rule numbers
function Base.parse(file, ::Type{RubiRules})
	regex = r"(Int.+) := (.+) \/; (.+)( \*\))?"
	# Capture the lines that fit the regex, skip the empty ones, put it all into
	# a big array. Gotta love Julia!
	caps = [i for i in match.(regex, readlines(file)) if !isnothing(i)]
	vlhs = lhs.(caps)
	vrhs = rhs.(caps)
	vgivens = givens.(caps)
	vcomments = iscommented.(caps)
	rubi = artifact"Rubi"
	path = relpath(file, rubi) # Fix this global?
	filename = splitpath(file)[end]
	# Assemble and return 📦 ➡
	[IntRuleCapture(pathname = path,
					filename = filename,
					rulenumber = 0,
					lhs = vlhs[i],
					rhs = vrhs[i],
					givens = vgivens[i],
					comment = vcomments[i]) for i in 1:length(caps)]
end

## Testing nothing's borked with finding the files and making a TreeFile of them
intfiles = rubidir["Rubi-4.16.1.0/Rubi/IntegrationRules"];
parsedstructs = FileTrees.load(intfiles) do file
	parse(string(path(file)), RubiRules)
end;
vstructs = reducevalues(vcat, parsedstructs);
@test 7032 == length(vstructs)

"""
Loads the RubiRules into a vector.
See `IntRuleCapture` for the fields.
We use FileTrees.jl, so starting Julia with multiple threads should speed this up,
but it feels like it could be faster ... ?
"""
function load(::Type{RubiRules})
	rubi = artifact"Rubi"
	intfiles = joinpath(rubi, "Rubi-4.16.1.0", "Rubi", "IntegrationRules")
	files = FileTree(intfiles)
	parsedstructs = FileTrees.load(files) do file
		parse(string(path(file)), RubiRules)
	end
	vstructs = reducevalues(vcat, parsedstructs)
	7032 == length(vstructs) || error("Please alert @miguelraz, something has 💥")
	vstructs
end

## Writing to JSON! At last!
# Defining this straight from the JSON3 documentation
StructTypes.StructType(::Type{IntRuleCapture}) = StructTypes.Struct();
# For debugging
#json = JSON3.write(vstructs);
"""
Write `vstructs` into a `intrules.json` file in Rubin.jl/src, where `vstructs`
is defined as
```julia
vstructs = load(RubiRules)
```
"""
function write(vstructs::Vector{IntRuleCapture}, ::Type{RubiRules})
	# hat tip to Jacob Quinn and the #data gang
	# Find the Rubin root dir in an OS independent way
	targetpath = joinpath(pkgdir(Rubin), "src", "rubirules.json")
	# JSON3 way to write to a file with 💃 pretty printing 💃
	open(targetpath, "w") do f
		JSON3.pretty(f, JSON3.write(vstructs))
		println(f)
	end
end

##################################################################
## Chapter 2:
#        The friggin' tests
#
# GOALS: ⚽
#1. Read MathematicaSyntaxTestSuite into a huge JSON Array
#2. Add metadata, facts, number, filename, test number
#3. Save to a single file.
#
## Note: They claim 72944 tests total in the website,
# but we must care about the ones that are commented
#
#                  💪 Currently at 72961 💪,
#                          with 1??        missing
#                            Wermer book missing
##################################################################

## Regex tests
# Tricky, gotta make sure we don't miss any of them!
inttests =
["{Sqrt[2*x + 1], x, 1, (1/3)*(1 + 2*x)^(3/2)}",
"(* {Sqrt[2*x + 1], x, 1, (1/3)*(1 + 2*x)^(3/2)} *)"];
inttestregex = r"{(.+),(.+),(.+),(.+)}";
res = match.(inttestregex, inttests);
@test res[1].captures[1] == "Sqrt[2*x + 1]"
@test res[1].captures[2] == " x"
@test res[1].captures[3] == " 1"
@test res[1].captures[4] == " (1/3)*(1 + 2*x)^(3/2)"
@test res[2].captures[1] == "Sqrt[2*x + 1]"
@test res[2].captures[2] == " x"
@test res[2].captures[3] == " 1"
@test res[2].captures[4] == " (1/3)*(1 + 2*x)^(3/2)"

## Header regex
# We'll likely need some metadata to alleviate pain when updating the codebase.
# This should help mitigate some of those aches.
headerregex = r"\(\*(.+)\*\)";
@test match(headerregex, "(*Integrands of the form x^m PolyLog[n, a x^q]*)").captures[1] == "Integrands of the form x^m PolyLog[n, a x^q]"

"""
 Sample input: {Sqrt[2*x + 1], x,   1, (1/3)*(1 + 2*x)^(3/2)}
               ^query		    ^var ^steps      ^optimal answer
# NOTE:
grep '(* {' | wc detects 284 commented cases at time of writing

# Fields:
	- filename
	- path
	- header
	- integrand
	- variable
	- steps
	- optimal
	- iscomment
"""
Base.@kwdef struct IntRuleTest
	filename::String =""
	path::String = ""
	integrand::String = ""
	variable::String = ""
	steps::String = "" # Yeah sure, this could be an Int but I'm 🏈 punting.
	optimal::String = ""
	iscomment::Bool = false
end

# header: string - to print testset and @async regions
# tests: Vec{IntRuleTest} - to write all the tests
Base.@kwdef struct IntRuleTestSection
	filename::String = ""
	path::String = ""
	header::String = ""
	tests::Vector{IntRuleTest} = IntRuleTest[]
end

"""
 1. Find first line that is not `skipme`, parse and push into the header
 2. Slurp the subsequent integrals, capture them into a vector of `IntRuleTest`s
 3. When you encounter a header that is different from the previous, put
		all the previous IntRuleTests into a `IntTestSectionCapture` array
	Headers: There is much metadata like `(* n > 0 *)` and `(*Integrands ... )`
	We don't wanna throw that away if possible
# TODO: finer metadata with the headers in the file
So that testset nesting can happen
"""
function Base.parse(file, ::Type{RubiTests})
	# TODO add nesting to subsubsections
	# This regex globs up everything (*InsideTheirCommentSyntax*)
	headerregex = r"\(\*(.+)\*\)"
	# This regex globs up the 4 entries in {a,format,like,this}
	# TODO Fix those 5 missing cases in issue
	# inttestregex = r"{(.+),(.+),(.+),(.+)?}"
	inttestregex = r"
    {(.+),                  # captures integrand
    \s?([a-z])\s?,          # captures variable (single letter)
    \s?((?:-?[0-9\s]+)|(?:If\[\$VersionNumber[>=<]{1,3}[-,\s\d]*\])) # captures steps (number or VersionNumber check)
    (?:,(.+))?}             # capturens optimal result if it exists
    "x
	rubi = artifact"RubiTests"
	path = relpath(file)

	header = ""
	tests = IntRuleTest[]
	for line in readlines(file)
		# Skip the `(* ::Subsection...` and `(* ::Subsubsection...` lines
		isempty(line) && continue
		startswith(line, "(* ::") && continue
		# If it's a test, push it
		if occursin(r"{", line) && occursin(r"}", line)
				# TODO Handle this better? We're losing information here
				# Perhaps adding a hierarchy to the section or subsection is good
			try
				m = match(inttestregex, line).captures
			catch
				@warn "inttestregex did not capture int test line"
				@warn line
				@warn file
				continue
			end
			m = match(inttestregex, line).captures

			# We need to know 7 things:
			# 1. integrand
			# 2. variable to integrate
			# 3. steps taken by rubi
			# 4. optimal answer
			# 5. if it's commented
			# 6. path
			# 7. filename
			integrand, variable, steps, optimal = (m...,)
			iscomment = startswith(line, "(* {")
			path = relpath(file, rubi) # Fix this global?
			filename = splitpath(file)[end]
			#@info "pushed to test"
			push!(tests, IntRuleTest(integrand = integrand, variable = variable,
									 path = path, filename =  filename,
									 steps = steps, optimal = string(optimal),
									 iscomment = iscomment))
		elseif startswith("(*", line) && endswith("*)", line)
		# 🆙 Update the header
			header = match(headerregex, line).captures[1]
		else
			continue
		end
	end
	# Remember! Push the last section!
	tests
end

## Test the test parsing function
# Because that is one chonky boii
#testfile = "/home/mrg/.julia/artifacts/1148cba18dae2f8939af8bc542233a48cc42cf19/MathematicaSyntaxTestSuite-4.16.0/5 Inverse trig functions/5.5 Inverse secant/5.5.2 Inverse secant functions.m";
testfile = "parsetest.m"
vtests = parse(testfile, RubiTests)
@test length(vtests) == 11
@test vtests[3].iscomment
@test vtests[9].iscomment

# TODO swap out the loads
function findtreefolder(::Type{RubiRules})
	rubi = artifact"Rubi"
	intfiles = joinpath(rubi, "Rubi-4.16.1.0", "Rubi", "IntegrationRules")
	files = FileTree(intfiles)
end
function findtreefolder(::Type{RubiTests})
	rubitests = artifact"RubiTests"
	intfiles = joinpath(rubitests, "MathematicaSyntaxTestSuite-4.16.0")
	files = FileTree(intfiles)
end
"""
	Loads the RubiRules into a vector.
See `IntTestSections` for the fields.
"""
function load(::Type{RubiTests})
	rubitests = artifact"RubiTests"
	intfiles = joinpath(rubitests, "MathematicaSyntaxTestSuite-4.16.0")
	fs = FileTree(intfiles)
	# Wester problems has an irregular format, see issue # 3
	rm(fs, r"LICENSE|README.md|Wester Problems.m")
	@info fs
	parsedstructs = FileTrees.load(fs) do file
		parse(string(path(file)), RubiTests)
	end
	vtests = reducevalues(vcat, parsedstructs)
	72961 == length(vtests) || error("Please alert @miguelraz, something has 💥")
	vtests
end

vtests = load(RubiTests);
@test length(vtests) == 72961

# For JSON3 writing
StructTypes.StructType(::Type{IntRuleTest}) = StructTypes.Struct();

"""
Write `vtests` into a `inttests.json` file in Rubin.jl/src, where `vtests`
is defined as
```julia
vtests = load(RubiRules)
```
"""
function write(vtests::Vector{IntRuleTest}, ::Type{RubiTests})
	# Hat tip to Jacob Quinn 🎩 and the #data gang
	# Find the Rubin root dir in an OS independent way
	targetpath = joinpath(pkgdir(Rubin), "src", "rubitests.json")
	# JSON3 way to write to a file with 💃 pretty printing 💃
	open(targetpath, "w") do f
		JSON3.pretty(f, JSON3.write(vtests))
		println(f)
	end
end

##################################################################
## Chapter 3:
#        Woflram -> Julia
#
# GOALS: ⚽
# 1. OUTPUT: tests/rubitests.jl with all the test files written in plain Julia
#	- Format shown below
# 2. OUTPUT: src/rubirules.jl with all the rules written in plain Julia
#   - Format: show below
#
# Approaches:
# 1. Make a small dict with Mathematica => Julia translations,
#	THEN search and replace every string on both sides
# 2. use a while loop with the replacings, if it doesn't change go on.
# OOOPS
#  - All fun and games until E^x and E^(x+2) showed up - unclear how to deal
#    with parens 💩
# 3. Use a library like RBNF.jl
# 4. DERP! See if someone else has solved it first!
# Sympy has a Mathematica parser - from there it should be much simple.
#  -
# Options:
#
# 1. Symbolics.jl vs Metatheory.jl .... vs Rust egraphs?
# # Consider how to propagate assumptions
# https://julialang.zulipchat.com/#narrow/stream/236639-symbolic-programming/topic/Integration/near/232350734
##################################################################

### Notes:
# Once I have a symbolic sympy expression, I can just turn it into a Julia
# Julia string with `str(y)`:
# sympy = pyimport("sympy")
# str = py"'Sin[x]^2'";
# mathparser = sympy.parsing.mathematica.mathematica
# y = maths(str)
# res = string(y) # 🚀
# @test "sin(x)^2" == res

## Parser hijinks:
# Ugh - the sympy.mathematica parser does not accept Int[ ... :(
# No worries 😅, we can pass in a PyDict with our prefs and have it work!
# function parsetosympy(s,d)
# 	s = pystring(PyObject(s))
#     f = sympy.parsing.mathematica.mathematica
#     f(s,d)
# end

## Parser hijinks: Sympy doesn't translate all functions :/
# Need to grep for the function names in rubirules.json and add them to the below dict
#mathfuncs = readlines("mathematicafunctions2.txt")
# mathfuncsdictionary = Dict(i * "*x]" => lowercase(i[1:end-1]) * "(*x)" for i in readlines("mathematicafunctions2.txt"))

# function mathfuncsdict()
# 	mathfuncs = readlines("mathematicafunctions2.txt")
# #	ps =
# end

# note the *x to slurp all args
sympydictfixes = PyDict(Dict("Int[*x]" => "integrate(*x)",
							 "ProductLog[*x]" => "productlog(*x)",
							 "Gamma[*x]" => "gamma(*x)",
							 # Following convention from SpecialFunctions.jl
							 "EllipticF[*x]" => "ellipf(*x)",
							 # TODO add Fresnel, FresnelS, FresnelC
							 #"E^x" => "exp(x)",
							 #"E^(*x)" => "exp(*x)",
							 "Abs[*x]" => "abs(*x)",
							 "AppellF1[*x]" => "appelf1(*x)",
							 "CannotIntegrate[*x]" => "cannotintegrate(*x)",
							 "CoshIntegral[*x]" => "coshintegral(*x)",
							 "CosIntegral[*x]" => "cosintegral(*x)",
							 "Defer[*x]" => "defer(*x)",
							 "Derivative[*x]" => "derivative(*x)",
							 "EllipticE[*x]" => "ellipe(*x)",
							 "EllipticF[*x]" => "ellipf(*x)",
							 "EllipticK[*x]" => "ellipk(*x)",
							 "EllipticPi[*x]" => "ellippi(*x)",
							 "Erf[*x]" => "erf(*x)",
							 "Erfc[*x]" => "erfc(*x)",
							 "Erfi[*x]" => "erfi(*x)",
							 "Expand[*x]" => "expand(*x)",
							 "ExpIntegralE[*x]" => "expintegrale(*x)",
							 "ExpIntegralEi[*x]" => "expintegralei(*x)",
							 #"f[*x]" => f(*x)",
							 #"F[*x]" => f(*x)",
							 "F0[*x]" => "f0(*x)",
							 "FresnelC[*x]" => "fresnelc(*x)",
							 "FresnelS[*x]" => "fresnels(*x)",
							 #"g[*x]" => ,
							 #"Gamma[*x]" => ,
							 "Hypergeometric2F1[*x]" => "hypergeometric2f1(*x)",
							 "HypergeometricPFQ[*x]" => "hypergeometricpfq(*x)",
							 #"If[*x]" => if,
							 #"Log[*x]" => "log,
							 "LogGamma[*x]" => "loggamma(*x)",
							 "LogIntegral[*x]" => "logintegral(*x)",
							 "PolyGamma[*x]" => "polygamma(*x)",
							 "PolyLog[*x]" => "polylog(*x)",
							 "ProductLog[*x]" => "productlog(*x)",
							 #"sdx[*x]" => ,
							 "SinhIntegral[*x]" => "sinhintegral(*x)",
							 "SinIntegral[*x]" => "sinintegral(*x)",
							 "Unintegrable[*x]" => "unintegrable(*x)",
							 "Zeta[*x]" => "zeta(*x)",
							 ));

funcdicts = Dict(
	"Abs[*x]" => "abs(*x)",
	"ActivateTrig[*x]" => "activatetrig(*x)",
	"AlgebraicFunctionQ[*x]" => "algebraicfunctionq(*x)",
	"Apart[*x]" => "apart(*x)",
	"AppellF1[*x]" => "appellf1(*x)",
	"ArcCos[*x]" => "arccos(*x)",
	"ArcCosh[*x]" => "arccosh(*x)",
	"ArcCot[*x]" => "arccot(*x)",
	"ArcCoth[*x]" => "arccoth(*x)",
	"ArcCsc[*x]" => "arccsc(*x)",
	"ArcCsch[*x]" => "arccsch(*x)",
	"ArcSec[*x]" => "arcsec(*x)",
	"ArcSech[*x]" => "arcsech(*x)",
	"ArcSin[*x]" => "arcsin(*x)",
	"ArcSinh[*x]" => "arcsinh(*x)",
	"ArcTan[*x]" => "arctan(*x)",
	"ArcTanh[*x]" => "arctanh(*x)",
	"AtomQ[*x]" => "atomq(*x)",
	"BesselJ[*x]" => "besselj(*x)",
	"BinomialDegree[*x]" => "binomialdegree(*x)",
	"BinomialMatchQ[*x]" => "binomialmatchq(*x)",
	"BinomialParts[*x]" => "binomialparts(*x)",
	"BinomialQ[*x]" => "binomialq(*x)",
	"Block[*x]" => "block(*x)",
	"CalculusFreeQ[*x]" => "calculusfreeq(*x)",
	"Cancel[*x]" => "cancel(*x)",
	"CannotIntegrate[*x]" => "cannotintegrate(*x)",
	"Coeff[*x]" => "coeff(*x)",
	"Coefficient[*x]" => "coefficient(*x)",
	"ComplexFreeQ[*x]" => "complexfreeq(*x)",
	"Complex[*x]" => "complex(*x)",
	"CosIntegral[*x]" => "cosintegral(*x)",
	"Cos[*x]" => "cos(*x)",
	"CoshIntegral[*x]" => "coshintegral(*x)",
	"Cosh[*x]" => "cosh(*x)",
	"Cot[*x]" => "cot(*x)",
	"Coth[*x]" => "coth(*x)",
	"Csc[*x]" => "csc(*x)",
	"Csch[*x]" => "csch(*x)",
	"D[*x]" => "d(*x)",
	"DeactivateTrig[*x]" => "deactivatetrig(*x)",
	"Denom[*x]" => "denom(*x)",
	"Denominator[*x]" => "denominator(*x)",
	"DerivativeDivides[*x]" => "derivativedivides(*x)",
	"Derivative[*x]" => "derivative(*x)",
	"Discriminant[*x]" => "discriminant(*x)",
	"Dist[*x]" => "dist(*x)",
	"Distrib[*x]" => "distrib(*x)",
	"DistributeDegree[*x]" => "distributedegree(*x)",
	"Divides[*x]" => "divides(*x)",
	"EllipticE[*x]" => "elliptice(*x)",
	"EllipticF[*x]" => "ellipticf(*x)",
	"EllipticPi[*x]" => "ellipticpi(*x)",
	"EqQ[*x]" => "eqq(*x)",
	"Erf[*x]" => "erf(*x)",
	"Erfc[*x]" => "erfc(*x)",
	"Erfi[*x]" => "erfi(*x)",
	"EulerIntegrandQ[*x]" => "eulerintegrandq(*x)",
	"EveryQ[*x]" => "everyq(*x)",
	"ExpIntegralE[*x]" => "expintegrale(*x)",
	"ExpIntegralEi[*x]" => "expintegralei(*x)",
	"Exp[*x]" => "exp(*x)",
	"ExpandExpression[*x]" => "expandexpression(*x)",
	"ExpandIntegrand[*x]" => "expandintegrand(*x)",
	"ExpandLinearProduct[*x]" => "expandlinearproduct(*x)",
	"ExpandToSum[*x]" => "expandtosum(*x)",
	"ExpandTrigExpand[*x]" => "expandtrigexpand(*x)",
	"ExpandTrigReduce[*x]" => "expandtrigreduce(*x)",
	"ExpandTrigToExp[*x]" => "expandtrigtoexp(*x)",
	"ExpandTrig[*x]" => "expandtrig(*x)",
	"Expand[*x]" => "expand(*x)",
	"Expon[*x]" => "expon(*x)",
	"Exponent[*x]" => "exponent(*x)",
	"F[*x]" => "f(*x)",
	"F_[*x]" => "f_(*x)",
	"Factor[*x]" => "factor(*x)",
	"FalseQ[*x]" => "falseq(*x)",
	"Floor[*x]" => "floor(*x)",
	"FracPart[*x]" => "fracpart(*x)",
	"FractionQ[*x]" => "fractionq(*x)",
	"FractionalPart[*x]" => "fractionalpart(*x)",
	"FreeFactors[*x]" => "freefactors(*x)",
	"FreeQ[*x]" => "freeq(*x)",
	"FresnelC[*x]" => "fresnelc(*x)",
	"FresnelS[*x]" => "fresnels(*x)",
	"FullSimplify[*x]" => "fullsimplify(*x)",
	"FunctionOfExponentialFunction[*x]" => "functionofexponentialfunction(*x)",
	"FunctionOfExponentialQ[*x]" => "functionofexponentialq(*x)",
	"FunctionOfExponential[*x]" => "functionofexponential(*x)",
	"FunctionOfLinear[*x]" => "functionoflinear(*x)",
	"FunctionOfLog[*x]" => "functionoflog(*x)",
	"FunctionOfQ[*x]" => "functionofq(*x)",
	"FunctionOfSquareRootOfQuadratic[*x]" => "functionofsquarerootofquadratic(*x)",
	"FunctionOfTrigOfLinearQ[*x]" => "functionoftrigoflinearq(*x)",
	"FunctionOfTrig[*x]" => "functionoftrig(*x)",
	"Function[*x]" => "function(*x)",
	"GCD[*x]" => "gcd(*x)",
	"G[*x]" => "g(*x)",
	"G_[*x]" => "g_(*x)",
	"Gamma[*x]" => "gamma(*x)",
	"GeQ[*x]" => "geq(*x)",
	"GeneralizedBinomialMatchQ[*x]" => "generalizedbinomialmatchq(*x)",
	"GeneralizedBinomialQ[*x]" => "generalizedbinomialq(*x)",
	"GeneralizedTrinomialDegree[*x]" => "generalizedtrinomialdegree(*x)",
	"GeneralizedTrinomialMatchQ[*x]" => "generalizedtrinomialmatchq(*x)",
	"GeneralizedTrinomialQ[*x]" => "generalizedtrinomialq(*x)",
	"GtQ[*x]" => "gtq(*x)",
	"H[*x]" => "h(*x)",
	"H_[*x]" => "h_(*x)",
	"HalfIntegerQ[*x]" => "halfintegerq(*x)",
	"Head[*x]" => "head(*x)",
	"Hold[*x]" => "hold(*x)",
	"HyperbolicQ[*x]" => "hyperbolicq(*x)",
	"Hypergeometric2F1[*x]" => "hypergeometric2f1(*x)",
	"HypergeometricPFQ[*x]" => "hypergeometricpfq(*x)",
	"IGeQ[*x]" => "igeq(*x)",
	"IGtQ[*x]" => "igtq(*x)",
	"ILeQ[*x]" => "ileq(*x)",
	"ILtQ[*x]" => "iltq(*x)",
	"Identity[*x]" => "identity(*x)",
	"If[*x]" => "if(*x)",
	"IndependentQ[*x]" => "independentq(*x)",
	"InertTrigFreeQ[*x]" => "inerttrigfreeq(*x)",
	"InertTrigQ[*x]" => "inerttrigq(*x)",
	"IntBinomialQ[*x]" => "intbinomialq(*x)",
	"IntHide[*x]" => "inthide(*x)",
	"IntLinearQ[*x]" => "intlinearq(*x)",
	"IntPart[*x]" => "intpart(*x)",
	"IntQuadraticQ[*x]" => "intquadraticq(*x)",
	"IntSum[*x]" => "intsum(*x)",
	"Int[*x]" => "int(*x)",
	"IntegerPart[*x]" => "integerpart(*x)",
	"IntegerQ[*x]" => "integerq(*x)",
	"IntegersQ[*x]" => "integersq(*x)",
	"IntegralFreeQ[*x]" => "integralfreeq(*x)",
	"Integral[*x]" => "integral(*x)",
	"Integrate[*x]" => "integrate(*x)",
	"InverseFunctionFreeQ[*x]" => "inversefunctionfreeq(*x)",
	"InverseFunctionOfLinear[*x]" => "inversefunctionoflinear(*x)",
	"InverseFunctionQ[*x]" => "inversefunctionq(*x)",
	"KnownCotangentIntegrandQ[*x]" => "knowncotangentintegrandq(*x)",
	"KnownSecantIntegrandQ[*x]" => "knownsecantintegrandq(*x)",
	"KnownSineIntegrandQ[*x]" => "knownsineintegrandq(*x)",
	"KnownTangentIntegrandQ[*x]" => "knowntangentintegrandq(*x)",
	"LeQ[*x]" => "leq(*x)",
	"LeafCount[*x]" => "leafcount(*x)",
	"LinearMatchQ[*x]" => "linearmatchq(*x)",
	"LinearPairQ[*x]" => "linearpairq(*x)",
	"LinearQ[*x]" => "linearq(*x)",
	"LogGamma[*x]" => "loggamma(*x)",
	"LogIntegral[*x]" => "logintegral(*x)",
	"Log[*x]" => "log(*x)",
	"LtQ[*x]" => "ltq(*x)",
	"MatchQ[*x]" => "matchq(*x)",
	"MemberQ[*x]" => "memberq(*x)",
	"MinimumMonomialExponent[*x]" => "minimummonomialexponent(*x)",
	"Mod[*x]" => "mod(*x)",
	"Module[*x]" => "module(*x)",
	"NeQ[*x]" => "neq(*x)",
	"NegQ[*x]" => "negq(*x)",
	"NiceSqrtQ[*x]" => "nicesqrtq(*x)",
	"NonfreeFactors[*x]" => "nonfreefactors(*x)",
	"NonsumQ[*x]" => "nonsumq(*x)",
	"NormalizeIntegrand[*x]" => "normalizeintegrand(*x)",
	"NormalizePowerOfLinear[*x]" => "normalizepoweroflinear(*x)",
	"NormalizePseudoBinomial[*x]" => "normalizepseudobinomial(*x)",
	"Not[*x]" => "not(*x)",
	"Numer[*x]" => "numer(*x)",
	"Numerator[*x]" => "numerator(*x)",
	"OddQ[*x]" => "oddq(*x)",
	"PerfectSquareQ[*x]" => "perfectsquareq(*x)",
	"PiecewiseLinearQ[*x]" => "piecewiselinearq(*x)",
	"PolyGCD[*x]" => "polygcd(*x)",
	"PolyGamma[*x]" => "polygamma(*x)",
	"PolyLog[*x]" => "polylog(*x)",
	"PolyQ[*x]" => "polyq(*x)",
	"PolynomialDivide[*x]" => "polynomialdivide(*x)",
	"PolynomialInQ[*x]" => "polynomialinq(*x)",
	"PolynomialInSubst[*x]" => "polynomialinsubst(*x)",
	"PolynomialQ[*x]" => "polynomialq(*x)",
	"PolynomialQuotient[*x]" => "polynomialquotient(*x)",
	"PolynomialRemainder[*x]" => "polynomialremainder(*x)",
	"PosQ[*x]" => "posq(*x)",
	"PowerOfLinearMatchQ[*x]" => "poweroflinearmatchq(*x)",
	"PowerOfLinearQ[*x]" => "poweroflinearq(*x)",
	"PowerQ[*x]" => "powerq(*x)",
	"PowerVariableExpn[*x]" => "powervariableexpn(*x)",
	"ProductLog[*x]" => "productlog(*x)",
	"ProductQ[*x]" => "productq(*x)",
	"PseudoBinomialPairQ[*x]" => "pseudobinomialpairq(*x)",
	"QuadraticMatchQ[*x]" => "quadraticmatchq(*x)",
	"QuadraticProductQ[*x]" => "quadraticproductq(*x)",
	"QuadraticQ[*x]" => "quadraticq(*x)",
	"QuotientOfLinearsParts[*x]" => "quotientoflinearsparts(*x)",
	"QuotientOfLinearsQ[*x]" => "quotientoflinearsq(*x)",
	"Quotient[*x]" => "quotient(*x)",
	"RationalFunctionExpand[*x]" => "rationalfunctionexpand(*x)",
	"RationalFunctionExponents[*x]" => "rationalfunctionexponents(*x)",
	"RationalFunctionQ[*x]" => "rationalfunctionq(*x)",
	"RationalQ[*x]" => "rationalq(*x)",
	"RemoveContent[*x]" => "removecontent(*x)",
	"ReplaceAll[*x]" => "replaceall(*x)",
	"Rt[*x]" => "rt(*x)",
	"Sec[*x]" => "sec(*x)",
	"Sech[*x]" => "sech(*x)",
	"ShowStep[*x]" => "showstep(*x)",
	"Sign[*x]" => "sign(*x)",
	"Simp[*x]" => "simp(*x)",
	"SimplerIntegrandQ[*x]" => "simplerintegrandq(*x)",
	"SimplerQ[*x]" => "simplerq(*x)",
	"SimplerSqrtQ[*x]" => "simplersqrtq(*x)",
	"SimplifyIntegrand[*x]" => "simplifyintegrand(*x)",
	"Simplify[*x]" => "simplify(*x)",
	"SinIntegral[*x]" => "sinintegral(*x)",
	"Sin[*x]" => "sin(*x)",
	"SinhIntegral[*x]" => "sinhintegral(*x)",
	"Sinh[*x]" => "sinh(*x)",
	"SplitProduct[*x]" => "splitproduct(*x)",
	"Sqrt[*x]" => "sqrt(*x)",
	"SubstForFractionalPowerOfLinear[*x]" => "substforfractionalpoweroflinear(*x)",
	"SubstForFractionalPowerOfQuotientOfLinears[*x]" => "substforfractionalpowerofquotientoflinears(*x)",
	"SubstForFractionalPowerQ[*x]" => "substforfractionalpowerq(*x)",
	"SubstForInverseFunction[*x]" => "substforinversefunction(*x)",
	"SubstFor[*x]" => "substfor(*x)",
	"Subst[*x]" => "subst(*x)",
	"SumQ[*x]" => "sumq(*x)",
	"SumSimplerQ[*x]" => "sumsimplerq(*x)",
	"Sum[*x]" => "sum(*x)",
	"Tan[*x]" => "tan(*x)",
	"Tanh[*x]" => "tanh(*x)",
	"Together[*x]" => "together(*x)",
	"TrigQ[*x]" => "trigq(*x)",
	"TrigSimplifyQ[*x]" => "trigsimplifyq(*x)",
	"TrigSimplify[*x]" => "trigsimplify(*x)",
	"TrinomialMatchQ[*x]" => "trinomialmatchq(*x)",
	"TrinomialQ[*x]" => "trinomialq(*x)",
	"TryPureTanSubst[*x]" => "trypuretansubst(*x)",
	"Unintegrable[*x]" => "unintegrable(*x)",
	"With[*x]" => "with(*x)",
	"Zeta[*x]" => "zeta(*x)",
	"cos[*x]" => "cos(*x)",
	"cot[*x]" => "cot(*x)",
	"csc[*x]" => "csc(*x)",
	"f[*x]" => "f(*x)",
	"f_[*x]" => "f_(*x)",
	"g[*x]" => "g(*x)",
	"g_[*x]" => "g_(*x)",
	"lst[*x]" => "lst(*x)",
	"sec[*x]" => "sec(*x)",
	"sin[*x]" => "sin(*x)",
	"tan[*x]" => "tan(*x)",
	"tmp[*x]" => "tmp(*x)",
	"trig_[*x]" => "trig_(*x)",
	"uu[*x]" => "uu(*x)",
);

# @test "integrate(x,x)" == parsetosympy("Int[x,x]",sympydictfixes)
# "integrate(x,x)"

## Parser hijinks:
# 1. Need to strip underscores and `x_Symbol`

### OOPS:
# It is clear that many of the "steps" field are being mangled horribly
# steps = [length(v.steps) != 2 for v in vtests]
# @test count(steps) == length(vtests)

## New parsing method ============================================================
# https://discourse.julialang.org/t/best-practices-for-converting-a-mathematica-expression-to-julia/18335
const sympy_parsing_mathematica = SymPy.PyCall.pyimport("sympy.parsing.mathematica")
mathematica2julia(s::AbstractString, substitutions::Pair{<:AbstractString,<:AbstractString}...) =
SymPy.walk_expression(sympy_parsing_mathematica."mathematica"(s, Dict(substitutions...)))

mathematica2julia(s::AbstractString) = mathematica2julia(s, m2j_subs...)

# Pull data from func_subs file to create substitution pairs
#= #TODO: functions to add
	- Derivative
	- Hypergeometric2F1
	- f'[x]
	- Double check functions that might come from SpecialFunctions.jl, etc.
=#
m2j_subs = [r[1] => r[2] for r in  eachrow(readdlm("func_subs.csv", ';', String))]

function m2j(s::AbstractString; verbose=true)
	try
		return mathematica2julia(s)
	catch e
		verbose && @warn "Couldn't parse:\n$(s)\n"
		# @error e
		return "Parse Error"
	end
end

@show m2j("Csc[x + Sin[x]]")
@show m2j("Sin[2*x]")
@show m2j("EllipticF[1, 0.5]")
@show m2j("Csc[x + Sin[x]]")
@show m2j("(2*I*Sqrt[2]*EllipticE[Pi/4 - (I*x)/2, 2] * Sqrt[Sinh[x]])/Sqrt[I*Sinh[x]]")
@show m2j("((-Gamma[-n, -Log[t]])*(-Log[t])^n)/Log[t]^n")
@show m2j("(c + d*x)^0*PolyGamma[n, a + b*x]")
@show m2j("(c + d*x)^1*ExpIntegralE[n, a + b*x]")
@show m2j("Hypergeometric2F1[1/2, 1/2, 3/2, 1/2)")
@show m2j("Derivative[n][f][x]")
@show m2j("f'[x]")

## Saving work as I go in a dataframe to make incremental changes
# Getting rid of any tests that don't work right now.
fname = "test/conversion_data.csv"
# df = DataFrame(vtests) # Initial creation
df = DataFrame(CSV.File(fname))
@show names(df)
@show size(df)

# filter!(x -> !x.iscomment, df);
# failstrgs = ["CannotIntegrate", "Unintegrable"]
# filter!(x -> !any(occursin.(failstrgs, x.optimal)), df);
# transform!(df, :integrand => ByRow(m2j) => :juliaintegrand) # slooooow
# @show size(df)
# CSV.write(fname, df)
# Save Point --------------------------------------------

# @show size(df)
# filter!(x -> x.juliaintegrand != "Parse Error", df);
# @show size(df)
# CSV.write(fname, df)
# Save Point --------------------------------------------

# filter!(x -> !occursin("Hypergeometric2F1", x.optimal), df); 	# Removes about 4000
# filter!(x -> !occursin("VersionNumber", x.optimal), df); 		# Removes about 200
# transform!(df, :optimal => ByRow(x -> m2j(x, verbose=false)) => :juliaoptimal) # slooooow
# filter!(x -> x.juliaoptimal != "Parse Error", df);				# Removes about 3500
# filter!(x -> x.juliaoptimal != "0", df);
# @show size(df)
# CSV.write(fname, df)
# Save Point --------------------------------------------

create_local_path(s) = joinpath("test", split(s, "\\")[2:end-1]...)
create_fname(s) = replace(split(s, "\\")[end], ".m"=>".jl")
transform!(df, :path => ByRow(create_local_path) => :outputpath)
transform!(df, :path => ByRow(create_fname) => :outputfname)

# Makes all the directories, only need to run once if the folders are deleted
# unique(df.outputpath) .|> mkpath

## Write a test files =======================================================

preamble = """
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions\n
"""

gdf = groupby(df, :outputfname);

function fix_syntax(s)
	s = replace(s, " "=> "") # Remove spaces
	s = replace(s, "__" => ".")
end

for ((fn,), group) in pairs(gdf)
	N,M = size(group)
    @info "$N\ttests\tCreating File: $fn"

	outpath = joinpath(first(group.outputpath), fn)

	open(outpath, "w") do io
		write(io, preamble)
		for t in eachrow(group)
			i = fix_syntax(t[:juliaintegrand])
			v = t[:variable]
			o = fix_syntax(t[:juliaoptimal])
			write(io, "@test integrate($i, $v) == :($o)\n")
		end
	end

end

## TODO make a smaller test Suite in a single file
# TODO: create a smarter way to created abbreviated tests than just first
outpath = joinpath("test", "abbreviated_tests.jl")
@info "Creating File: $outpath"

N = 10
open(outpath, "w") do io
	write(io, preamble)
	for ((fn,), group) in pairs(gdf)
		setname = first(group.outputfname)
		@info "Creating test set: $setname"
		# Begin test set
		write(io, "@testset \"$setname\" begin\n")

		for t in eachrow(first(group,N))
			i = fix_syntax(t[:juliaintegrand])
			v = t[:variable]
			o = fix_syntax(t[:juliaoptimal])
			write(io, "\t@test integrate($i, $v) == :($o)\n")
		end

		# End of test set
		write(io, "end\n\n")
	end
end

##################################################################
## Chapter 4: Parsing Utility Functions
#
# Mathematica pattern strings:
# 	⋅ _			any single expression
# 	⋅ x_		any single expression, to be named x
# 	⋅ __		any sequence of one or more expressions
# 	⋅ x__		sequence named x
# 	⋅ x__h		sequence of expressions, all of whose heads are h
# 	⋅ ___		any sequence of zero or more expressions
# 	⋅ x___		sequence of zero or more expressions named x
# 	⋅ x___h		sequence of zero or more expressions, all of whose heads are h
#
# Mathematica heads:
# 	⋅ x_h		an expression with head h
# 	⋅ x_Integer	an integer
# 	⋅ x_Real	an approximate real number
# 	⋅ x_Complex	a complex number
# 	⋅ x_List	a list
# 	⋅ x_Symbol	a symbol
#
# Conditions:
# 	⋅ patt/;test		is a pattern which matches only if the evaluation of test yields True.
# 	⋅ lhs:>rhs/;test	represents a rule which applies only if the evaluation of test yields True.
# 	⋅ lhs:=rhs/;test	is a definition to be used only if test yields True.
##################################################################
struct RubiUtility <: AbstractRubiParser end

"""
 Sample input: IntHide[u_,x_Symbol] := Block[{\$ShowSteps=False,\$StepCounter=Null}, Int[u,x]]

# Fields:
	- name
	- args
	- eqtype
	- defintion
"""
Base.@kwdef struct UtilFuncMethod
	name::String = ""
	args::String = ""
	eqtype::String = ""
	defintion::String = ""
end

# Create util function type from regex matches automatically
UtilFuncMethod(re::RegexMatch) = UtilFuncMethod(re.captures...)
Base.show(io::IO, m::UtilFuncMethod) = print(io, "$(m.name)$(m.args)$(m.eqtype)$(m.defintion)")

# Needs for JSON3 to write a file
StructTypes.StructType(::Type{UtilFuncMethod}) = StructTypes.StringType()

"""
Sample input:
	(* ::Subsection::Closed:: *)
	(*TrigQ[u]*)

	TrigQ::usage = "If u is an expression of the form F[v] where F is a circular trig function, TrigQ[u] returns True; else it returns False.";
	\$TrigFunctions = {Sin, Cos, Tan, Cot, Sec, Csc};
	TrigQ[u_] := MemberQ[\$TrigFunctions, If[AtomQ[u],u,Head[u]]]

# Fields:
	- funcname
	- funlist
	- usage
	- methods
"""
Base.@kwdef struct IntUtilityFunc
	funcname::String = ""
	funlist::String = ""
	usage::String = ""
	methods::Vector{UtilFuncMethod} = []
	section::String = ""
end

# Pretty printing of functions with methods
function Base.show(io::IO, f::IntUtilityFunc)
	println(io, "$(f.funcname)")
	length(f.usage) > 0 && println("\tUsage:\t$(f.usage)")
	# f.methods |> enumerate .|> t -> println("  Method $(t[1]):  $(t[2])")
end

StructTypes.StructType(::Type{IntUtilityFunc}) = StructTypes.Struct()

function write(utilfuncs::Vector{IntUtilityFunc}, ::Type{RubiUtility})
	targetpath = joinpath(pkgdir(Rubin), "src", "rubiutils.json")
	open(targetpath, "w") do f
		JSON3.pretty(f, JSON3.write(utilfuncs))
		println(f)
	end
end

##
# TODO: parse by section and subsection???
# 	Requires splitting methods that doesn't eat the secregex
# 	Handeling Item sections?
#
# TODO: Only 4 functions fail to match
# 	⋅ IntegerPowerQ[u]					- no brackets
# 	⋅ FractionalPowerQ[u]				- no brackets
# 	⋅ ExpandBinomial[u,x]				- broken up into sections???
# 	⋅ DeactivateInverseHyperbolic[u,x] 	- commented out
#
function Base.parse(file, ::Type{RubiUtility})
	secregex = r"\(\* ::[S].*\*\)" # match get eaten during the split operation

	usageregex = r"""\S+::usage[\s="]+(.*)";"""
	# Matches the usage lines and ...   ^ captures usage string

	funcheaderregex = r"\(\*(\w*\[.*\])\*\)\s?"
	# Matches function header & ^ captures function name

	# funlistregex = r"(^\$\w+)(\s?=\s?)(\{.*\};)"
	funlistregex = r"\$\w+\s?=\s?\{(.*)\};" # w/o captures
	# Matches lists of functions

	# sec1nameregex = r"\(\* ::Section.*\n\(\*([\w\s]+)\*\)"
	# sec2nameregex = r"\(\* ::Subsection.*\n\(\*([\w\s]+)\*\)"
	secnameregex = r"\(\*([\w\s]{2,})\*\)"

	# funclineregex = r"\n([A-Za-z0-9]+)(\[[:A-Za-z0-9_,\^\*\+\[\]]*\])(\s*[:=]+\s*)((?:.+\n)+)" # oneliner
	funclineregex = r"
	\n([A-Za-z0-9]+)					# captures method name staring on new line
	(\[[:A-Za-z0-9_,\^\*\+\[\]]*\])		# captures methods arguments
	(\s*[:=]+\s*)						# captures equality operator (:= or =)
	((?:.+\n)+)							# captures multile lines of definition
	"x

	utilcontents = read(file, String);
	sections = split(utilcontents, secregex, keepempty=false)

	secname = ""
	utilfuncs = IntUtilityFunc[]
	for sec in sections[2:end]
		snm = match(secnameregex, sec) 		# Update section name
		if !isnothing(snm)
			secname = snm.captures[1]
		end

		fhm = match(funcheaderregex, sec)	# Function header with name
		isnothing(fhm) && continue
		funcname = fhm.captures[1]

		flm = match(funlistregex, sec)		# Function list
		funlist = isnothing(flm) ? "" : string(flm.match)

		usm = match(usageregex, sec)		# Usage string
		usage = usm isa Nothing ? "" : usm.captures[1]

		# Get function methods
		!occursin(funclineregex, sec) && (@warn "No methods detected for: $funcname"; continue)
		# create a list of methods using dispatch on Regex
		methods = UtilFuncMethod.(eachmatch(funclineregex, sec))

		# Push the utility function to the list
		push!( utilfuncs, IntUtilityFunc(
				funcname,
				funlist,
				usage,
				methods,
				secname
			)
		)
	end
	utilfuncs
end

function load(::Type{RubiUtility})
	rubi = artifact"Rubi"
	utilityfile = joinpath(rubi, "Rubi-4.16.1.0", "Rubi", "IntegrationUtilityFunctions.m")
	parse(utilityfile, RubiUtility)
end

utilfuncs = load(RubiUtility);
@test length(utilfuncs) >= 92

# Writing the JSON file
# json = write(utilfuncs, RubiUtility);


text:
function jacobian_transpose_ik!(state::MechanismState,
                               body::RigidBody,
                               point::Point3D,
                               desired::Point3D;
                               α = 0.1,
                               iterations = 100,
                               visualize = false,
                               mvis = "")
    mechanism = state.mechanism
    world = root_frame(mechanism)

    # Compute the joint path from world to our target body
    p = path(mechanism, root_body(mechanism), body)
    # Allocate the point jacobian (we'll update this in-place later)
    Jp = point_jacobian(state, p, transform(state, point, world))

    q = copy(configuration(state))

    for i in 1:iterations
        # Update the position of the point
        point_in_world = transform(state, point, world)
        # Update the point's jacobian
        point_jacobian!(Jp, state, p, point_in_world)
        # Compute an update in joint coordinates using the jacobian transpose
        Δq = α * Array(Jp)' * (transform(state, desired, world) - point_in_world).v
        # Apply the update
        q .= configuration(state) .+ Δq
        set_configuration!(state, q)
        visualize && set_configuration!(mvis, q)
    end
    q
end

ext: jl

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
include("classic.jl")
using Groebner
using Test
using AbstractAlgebra

using Logging
Logging.global_logger(Logging.ConsoleLogger(Logging.Warn))

###############################################################################

zerobuch()

fs = Groebner.rootn(3)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
@test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_NF BUCH_DISCARDED BUCH_REDUCED

###############################################################################

zerobuch()

fs = Groebner.rootn(4)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
@test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_NF BUCH_DISCARDED BUCH_REDUCED

###############################################################################

zerobuch()

fs = Groebner.change_ordering(Groebner.noonn(3), :degrevlex)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
@test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_NF BUCH_DISCARDED BUCH_REDUCED

###############################################################################

zerobuch()

fs = Groebner.change_ordering(Groebner.noonn(4, ground=GF(2^31-1)), :degrevlex)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
@test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_NF BUCH_DISCARDED BUCH_REDUCED

###############################################################################

zerobuch()

fs = Groebner.change_ordering(Groebner.kinema(ground=GF(2^31-1)), :degrevlex)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
@test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_SPOLYS BUCH_DISCARDED BUCH_SPOLYS-BUCH_DISCARDED

###############################################################################

zerobuch()

fs = Groebner.change_ordering(Groebner.sparse5(), :degrevlex)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
# @test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_NF BUCH_DISCARDED BUCH_REDUCED

###############################################################################

zerobuch()

fs = Groebner.change_ordering(Groebner.eco5(), :degrevlex)
gb = Groebner.groebner(fs)

sign_gb = groebner_basis(fs)

println(sign_gb)
# @test Groebner.isgroebner(sign_gb)
@test Groebner.groebner(sign_gb) == gb
@error "" BUCH_NF BUCH_DISCARDED BUCH_REDUCED


text:
using Plots
using ForwardDiff
using Roots
include("helper_functions.jl")

begin
    λ1 = 0.5
    λ2 = 1
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    q = 0.5
    qs = 0:0.001:1
end

# Test 1: Do the derivatives agree?
# function adcost(p, q, λ1, λ2, α, μ, ν)
#     (λ1 * p) / (λ1 * p - ν(α) + λ2 * q)^2 -
#     1/(μ(α) - λ1 + λ1 * p) -
#     1/(λ1 * p - ν(α) + λ2 * q) +
#     (λ1 * (p - 1))/(μ(α) - λ1 + λ1 * p)^2
# end

# dcost = ∂cost_∂p

ps = 0:0.001:1
plot(ps, ∂cost_∂p.(ps, q, parameters...))
plot!(ps, diff_cost.(ps, q, parameters...))
plot!(ylims = (-500, 500))

# @benchmark adcost(0.5, q, parameters...)
# @benchmark dcost(0.5, q, parameters...)

# Test 2: Do the points where derivative is zero agree? (aka pstar)
sol, _ = first.(find_valid_roots(p -> ∂cost_∂p(p, q, parameters...)))
# function apstar(q, λ1, λ2, α, μ, ν)
#     p = (λ1*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) - μ(α)*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) - ν(α)*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) - λ1*ν(α) + 2*μ(α)*ν(α) + λ2*q*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) + λ1*λ2*q - 2*λ2*μ(α)*q)/(λ1*μ(α) - λ1*ν(α) + λ1*λ2*q) |> real
#     return clamp(p, 0, 1)
# end
asol = pstar(q, parameters...)

diff_cost(sol, q, parameters...)
diff_cost(asol, q, parameters...)

# Test 3: Do the pstars agree?
# function pstar(q, λ1, λ2, α, μ, ν)
#     roots = find_valid_roots(p -> ∂cost_∂p(p, q, λ1, λ2, α, μ, ν))
#     if length(roots[1]) > 1
#         root_values, root_is_minimum =
#             classify_roots(p -> ∂cost_∂p(p, q, λ1, λ2, α, μ, ν), roots)
#         return (root_values[root_is_minimum])[1]
#     else
#         return (roots[1])[1]
#     end
# end

qs = 0:0.0001:1
ps = numerical_pstar.(qs, parameters...)
a_ps = pstar.(qs, parameters...)

response_plot = plot(qs, ps, label = "autodiff")
plot!(qs, a_ps, label = "analytic")
plot_assumptions!(response_plot, [ps, a_ps], qs, parameters...)
plot!(response_plot, legend = :outertopleft, size = (690, 400))

costs = cost.(ps, qs, parameters...)
a_costs = cost.(a_ps, qs, parameters...)

cost_plot = plot(qs, costs, label = "autodiff")
plot!(qs, a_costs, label = "analytic")
plot!(ylims = (-50, 50))
plot_assumptions!(cost_plot, [ps, a_ps], qs, parameters..., height = 10)

traffic_plot = plot(
    title = "Traffic at each server for varying q",
    xlabel = "Opponent probability (q)",
    ylabel = "Intensity of arrivals",
    legend = :topleft,
    xlims = (0, 1),
    ylims = ylims = (0, Inf)
    )

plot!(qs, λ1 .* a_ps .+ λ2 .* qs, label = "ν(α)")
hline!([ν.(α)], label = "ν(α) limit", ls = :dash, lw = 1.5)

plot!(qs, λ1 .* (1 .- a_ps), label = "μ(α)")
hline!([μ(α)], label = "μ(α) limit", ls = :dash, lw = 1.5)

intersection_index = last(findmax(λ1 .* a_ps .+ λ2 .* qs .> ν(α)))
x_intersection = qs[intersection_index]
y_intersection_max = maximum([
    (λ1 .* (1 .- a_ps))[intersection_index],
    (λ1 .* a_ps .+ λ2 .* qs)[intersection_index]]
    )
plot!([x_intersection, x_intersection], [0, y_intersection_max], lw = 1, lc = :maroon, label = "")

traffic_plot_anim = @animate for α ∈ vcat(collect(0.5:0.005:2.1), collect(reverse(0.5:0.005:2.1)))
    a_ps = apstar.(qs, (λ1, λ2, α, μ, ν)...)
    traffic_plot = plot(
    title = "Traffic at each server for varying qs at α = $(dp_string(α, 3))",
    # xlabel = "Opponent probability (q)",
    ylabel = "Intensity of arrivals",
    legend = :topleft,
    xlims = (0, 1),
    ylims = ylims = (0, Inf),
    xticks = 0:0.2:1,
    )

    plot!(qs, λ1 .* a_ps .+ λ2 .* qs, label = "ν(α)")
    hline!([ν.(α)], label = "ν(α) limit", ls = :dash, lw = 1.5)

    plot!(qs, λ1 .* (1 .- a_ps), label = "μ(α)")
    hline!([μ(α)], label = "μ(α) limit", ls = :dash, lw = 1.5)

    does_intersect, intersection_index = findmax(λ1 .* a_ps .+ λ2 .* qs .> ν(α))
    x_intersection = qs[intersection_index]
    y_intersection_max = maximum([
        (λ1 .* (1 .- a_ps))[intersection_index],
        (λ1 .* a_ps .+ λ2 .* qs)[intersection_index]]
        )
    plot!([x_intersection, x_intersection], [0, y_intersection_max], ls = :dash, lw = 1, lc = :black, label = "")

    pstar_plot = plot(
        qs,
        a_ps,
        lc = :black,
        label = "Best Response (p⋆)",
        xlabel = "Opponent probability (q)",
        xticks = 0:0.2:1,
        ylabel = "Response probability (p)",
        xlims = (0, 1),
        ylims = (0, 1))
    plot!(
        [x_intersection, x_intersection],
        [0, a_ps[intersection_index]],
        label = "",
        ls = :dash,
        lc = :black)

    if does_intersect
        plot!(
            qs[1:intersection_index],
            zeros(length(qs[1:intersection_index])),
            fillrange = a_ps[1:intersection_index],
            fc = :green,
            fa = 0.35,
            label = "BR Feasibile")
        plot!(
            qs[intersection_index + 1:end],
            zeros(length(qs[intersection_index + 1:end])),
            fillrange = a_ps[intersection_index + 1:end],
            fc = :red,
            fa = 0.35,
            label = "BR Infeasibile")
    else
        plot!(
            qs,
            zeros(length(qs)),
            fillrange = a_ps,
            fc = :green,
            fa = 0.35,
            label = "BR Feasibile")
        # Have to put something here so it doesn't change the legend
        plot!(
            [1, 2],
            zeros(2),
            fillrange = [2, 3],
            fc = :red,
            fa = 0.35,
            label = "BR Infeasibile")
    end

    plot(traffic_plot, pstar_plot, layout = (2, 1), size = (800, 800))
end

gif(traffic_plot_anim, "traffic_plot_anim.gif", fps = 30)


text:
export
    WignerFunction,
    WignerSurface

#=
    Wigner function by Laguerre Polynominal
=#

abstract type CreateWignerMethod end
struct Load𝐖 <: CreateWignerMethod end
struct Calc𝐖 <: CreateWignerMethod end

function wigner(m::Integer, n::Integer, x::Vector{<:Real}, p::Vector{<:Real})
    w = gaussian_function(x, p) .*
        coefficient_of_wave_function(m, n) .*
        z_to_power(m, n, x, p) .*
        laguerre(m, n, x, p)

    return w
end

function create_wigner(
    m_dim::Integer,
    n_dim::Integer,
    x_range::AbstractRange,
    p_range::AbstractRange,
    ::Type{Calc𝐖}
)
    𝐰 = Array{ComplexF64,4}(undef, m_dim, n_dim, length(x_range), length(p_range))
    @sync for m in 1:m_dim
        Threads.@spawn for n in 1:n_dim
            𝐰[m, n, :, :] = wigner(m, n, collect(x_range), collect(p_range))
        end
    end

    return 𝐰
end

function create_wigner(
    m_dim::Integer,
    n_dim::Integer,
    x_range::AbstractRange,
    p_range::AbstractRange,
    bin_path::String,
    ::Type{Load𝐖}
)
    return load_𝐰(m_dim, n_dim, x_range, p_range, bin_path)
end

function create_wigner(
    m_dim::Integer,
    n_dim::Integer,
    x_range::AbstractRange,
    p_range::AbstractRange,
)
    bin_path = gen_wigner_bin_path(m_dim, n_dim, x_range, p_range)

    if isfile(bin_path)
        return create_wigner(m_dim, n_dim, x_range, p_range, bin_path, Load𝐖)
    end

    𝐰 = create_wigner(m_dim, n_dim, x_range, p_range, Calc𝐖)
    save_𝐰(bin_path, 𝐰)

    return 𝐰
end

mutable struct WignerFunction{T<:Integer, U<:AbstractRange}
    m_dim::T
    n_dim::T
    x_range::U
    p_range::U
    𝐰::Array{ComplexF64,4}

    function WignerFunction(
        m_dim::T,
        n_dim::T,
        x_range::U,
        p_range::U
    ) where {T<:Integer, U<:AbstractRange}
        !check_argv(m_dim, n_dim, x_range, p_range) && throw(ArgumentError)

        𝐰 = create_wigner(m_dim, n_dim, x_range, p_range)
        return new{T, U}(m_dim, n_dim, x_range, p_range, 𝐰)
    end
end

"""
    WignerFunction(x_range::AbstractRange, p_range::AbstractRange; dim=DIM)

Declare a generalized Wigner function with `x` range and `p` range

* `x_range`: range of position.
* `p_range`: range of momentum.
* `dim`: Maximum photon number for truncate, default is $DIM.

# Examples
```julia-repl
julia> state = VacuumState();

julia> wf = WignerFunction(-10:0.1:10, -10:0.1:10);

julia> w = wf(state);
```
"""
function WignerFunction(x_range::AbstractRange, p_range::AbstractRange; dim=DIM)
    return WignerFunction(dim, dim, x_range, p_range)
end

struct WignerSurface{T<:AbstractRange}
    x_range::T
    p_range::T
    𝐰_surface::Matrix{Float64}

    function WignerSurface(
        x_range::T,
        p_range::T,
        𝐰_surface::Matrix{Float64}
    ) where {T<:AbstractRange}
        return new{T}(x_range, p_range, 𝐰_surface)
    end
end

function (wf::WignerFunction)(state::StateMatrix)
    𝛒 = state.𝛒

    𝐰_surface = Matrix{Float64}(undef, length(wf.x_range), length(wf.p_range))
    @sync for i in 1:length(wf.x_range)
        Threads.@spawn for j in 1:length(wf.p_range)
            𝐰_surface[i, j] = real(sum(𝛒 .* wf.𝐰[:, :, i, j]))
        end
    end

    return WignerSurface(wf.x_range, wf.p_range, 𝐰_surface)
end

(wf::WignerFunction)(state::StateVector) = wf(StateMatrix(state))


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Compat, Random
include("convexpenalties.jl")

results_template = DataFrame(n=Int[],p=Int[], m=Int[], k=Int[], theSeed=Int[], gamma=Real[], lambda_rnk=Real[], lambda_nn=Real[],
gap_persp=Real[], gap_dcl=Real[], t_persp=Real[], t_dcl=Real[], t_nn=Real[], acc_persp=Real[], acc_dcl=Real[], acc_nn=Real[],
rnk_persp=Int[], rnk_dcl=Int[], rnk_nn=Int[])

results_exact=similar(results_template, 0)

for ARG in ARGS

array_num = parse(Int, ARG)
seeds=1:20

ns=10:10:200
k=10
gamma=1e6

lambda_rnk=1.0
lambda_nn=1.0

# n_obs_oos=1000
theSeed=seeds[(array_num-1)÷20+1]
n=ns[(array_num-1)%20+1]
m=n
p=n
Random.seed!(theSeed)
X=rand(m,p)

U_true=rand(p,k)
V_true=rand(n,k)
beta_true=U_true*V_true'
E=0.05*rand(m,n)

Y=X*beta_true.+E

for t in 1:1

        t_persp=@elapsed beta_est,ofv_dual, ofv_primal,=solvepersprelax(X, Y, m, n, p, 0.279883043, gamma)
        @show t_persp
        @show gap_persp=100*abs.(ofv_primal-ofv_dual)/ofv_primal
        @show rel_error_persp=norm(beta_est.-beta_true)/norm(beta_true)
        @show ofv_primal
        @show 1.0/(2.0*gamma)*norm(beta_est)^2

        u,v,=svd(beta_est)
        @show rnk_persp=sum(v.>1e-4)

        t_dcl=@elapsed beta_est,ofv_dual, ofv_primal,= solvepersprelax_DCL(X, Y, m, n, p, 0.01215675, gamma)
        @show gap_dcl=100*abs.(ofv_primal-ofv_dual)/ofv_primal

        @show rel_error_dcl=norm(beta_est.-beta_true)/norm(beta_true)

        u,v,=svd(beta_est)
        @show rnk_dcl=sum(v.>1e-4)

        t_nn=@elapsed beta_est,=solveNNrelax(X, Y, m, n, p, 0.000573504, gamma)

        @show rel_error_nn=norm(beta_est.-beta_true)/norm(beta_true)
        u,v,=svd(beta_est)
        @show rnk_nn=sum(v.>1e-4)

        push!(results_exact, [n, p, m, k, theSeed, gamma, lambda_rnk, lambda_nn,
        gap_persp, gap_dcl, t_persp, t_dcl, t_nn, rel_error_persp, rel_error_dcl, rel_error_nn,
        rnk_persp, rnk_dcl, rnk_nn])

end
CSV.write("createFig2_times.csv", results_exact, append=true)

end

max_stars_repo_path: Matrix Regression/createFigure2.jl

max_issues_repo_path: Matrix Regression/createFigure2.jl

max_forks_repo_path: Matrix Regression/createFigure2.jl


text:
using Colors

function referencespacecache(reference_space_key="annotation/ccf_2017"; resolution=25, manifest=referencespacemanifest)
    reference_space_cache.ReferenceSpaceCache(resolution, reference_space_key, manifest=manifest)
end

function getstructuretree(space=referencespacecache(), id=1) # Default id is mouse brain atlas
    id = Int(id)
    space.get_structure_tree(structure_graph_id=id)
end

function getstructurecolor(id::Union{Number, Missing}; tree = getstructuretree())
    !ismissing(id) || return RGB(0.0, 0.0, 0.0)
    id = Int(id)
    d = tree.get_structures_by_id([id])[1]
    !isnothing(d) || return RGB(0.0, 0.0, 0.0)
    c = RGB((d["rgb_triplet"]./256)...)
end
function getstructurecolors(ids::AbstractVector)
    tree = getstructuretree()
    getstructurecolor.(ids; tree)
end

function getstructuretreedepth(id) # How many parents does this structure have (+1)?
    id = Int(id)
    tree = getstructuretree()
    d = tree.get_structures_by_id([id])[1]
    length(d["structure_id_path"])
end

function getstructurename(id::Number)
    id = Int(id)
    tree = getstructuretree()
    d = tree.get_structures_by_id([id])[1]
    d["name"]
end
function getstructurename(acronym::String)
    tree = getstructuretree()
    d = tree.get_structures_by_acronym([acronym])[1]
    d["name"]
end

function getstructureid(acronym::String)
    tree = getstructuretree()
    d = tree.get_structures_by_acronym([acronym])[1]
    d["id"]
end

function getallstructureids(args...)
    t = getstructuretree(args...)
    d = t.get_name_map()
    ids = keys(d)
end

#! Can do the rest of this dict by eval?

function buildreferencespace(tree=getstructuretree(), annotation=getannotationvolume(), resolution=(25,25,25))
    if annotation isa Tuple
        annotation = annotation[1]
    end
    rsp = reference_space.ReferenceSpace(tree, annotation, resolution)
end

function buildstructuremask(id, referencespace=buildreferencespace())
    referencespace.make_structure_mask([id])
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

uniform_cdf(x::Real, a::Real, b::Real) = (x - a) / (b - a)
uniform_invcdf(u::Real, a::Real, b::Real) = (b - a) * u + a
uniform_invcdf_ladj(u::Real, a::Real, b::Real) = log(b - a)

# ToDo: Cache inv(b - a)
struct UniformCDFTrafo{T<:Real} <: VariateTransform{Univariate,UnitSpace,MixedSpace}
    a::T
    b::T
end

ValueShapes.varshape(trafo::UniformCDFTrafo) = ScalarShape{Real}()

target_space(trafo::UniformCDFTrafo) = UnitSpace()
source_space(trafo::UniformCDFTrafo) = MixedSpace()

function apply_vartrafo_impl(trafo::UniformCDFTrafo, src_v::Real, prev_ladj::Real)
    trg_v = uniform_cdf(src_v, trafo.a, trafo.b)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = - uniform_invcdf_ladj(trg_v, trafo.a, trafo.b)
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

function apply_vartrafo_impl(inv_trafo::InverseVT{<:UniformCDFTrafo}, src_v::Real, prev_ladj::Real)
    trafo = inv_trafo.orig
    trg_v = uniform_invcdf(src_v, trafo.a, trafo.b)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = uniform_invcdf_ladj(src_v, trafo.a, trafo.b)
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

logistic_cdf(x::Real, mu::Real, theta::Real) = inv(exp(-(x - mu) / theta) + one(x))
logistic_invcdf(u::Real, mu::Real, theta::Real) = log(u / (one(u) - u)) * theta + mu
logistic_invcdf_ladj(u::Real, mu::Real, theta::Real) = -log((u - u^2) / theta)

# ToDo: Cache inv(theta)
struct LogisticCDFTrafo{T<:Real} <: VariateTransform{Univariate,UnitSpace,InfiniteSpace}
    mu::T
    theta::T
end

ValueShapes.varshape(trafo::LogisticCDFTrafo) = ScalarShape{Real}()

target_space(trafo::LogisticCDFTrafo) = UnitSpace()
source_space(trafo::LogisticCDFTrafo) = InfiniteSpace()

function apply_vartrafo_impl(trafo::LogisticCDFTrafo, src_v::Real, prev_ladj::Real)
    trg_v = logistic_cdf(src_v, trafo.mu, trafo.theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = - logistic_invcdf_ladj(trg_v, trafo.mu, trafo.theta)
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

function apply_vartrafo_impl(inv_trafo::InverseVT{<:LogisticCDFTrafo}, src_v::Real, prev_ladj::Real)
    trafo = inv_trafo.orig
    trg_v = logistic_invcdf(src_v, trafo.mu, trafo.theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = logistic_invcdf_ladj(src_v, trafo.mu, trafo.theta)
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

exponential_cdf(x::Real, theta::Real) = one(x) - exp(x / -theta)
exponential_invcdf(u::Real, theta::Real) = -theta * log(one(u) - u)
exponential_cdf_ladj(x::Real, theta::Real) = (x / -theta) - log(theta)

# ToDo: Cache inv(theta) and log(theta)
struct ExponentialCDFTrafo{T<:Real} <: VariateTransform{Univariate,UnitSpace,MixedSpace}
    theta::T
end

ValueShapes.varshape(trafo::ExponentialCDFTrafo) = ScalarShape{Real}()

target_space(trafo::ExponentialCDFTrafo) = UnitSpace()
source_space(trafo::ExponentialCDFTrafo) = MixedSpace()

function apply_vartrafo_impl(trafo::ExponentialCDFTrafo, src_v::Real, prev_ladj::Real)
    trg_v = exponential_cdf(src_v, trafo.theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = exponential_cdf_ladj(src_v, trafo.theta)
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

function apply_vartrafo_impl(inv_trafo::InverseVT{<:ExponentialCDFTrafo}, src_v::Real, prev_ladj::Real)
    trafo = inv_trafo.orig
    trg_v = exponential_invcdf(src_v, trafo.theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = - exponential_cdf_ladj(trg_v, trafo.theta)
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

# ToDo: Cache inv(theta)
struct ScaledLogTrafo{T<:Real} <: VariateTransform{Univariate,InfiniteSpace,MixedSpace}
    theta::T
end

ValueShapes.varshape(trafo::ScaledLogTrafo) = ScalarShape{Real}()

target_space(trafo::ScaledLogTrafo) = InfiniteSpace()
source_space(trafo::ScaledLogTrafo) = MixedSpace()

function apply_vartrafo_impl(trafo::ScaledLogTrafo, src_v::Real, prev_ladj::Real)
    trg_v = log(src_v / trafo.theta)
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = - abs(trg_v + log(trafo.theta))
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

function apply_vartrafo_impl(inv_trafo::InverseVT{<:ScaledLogTrafo}, src_v::Real, prev_ladj::Real)
    trafo = inv_trafo.orig
    trg_v = exp(src_v) * trafo.theta
    if isnan(prev_ladj)
        var_trafo_result(trg_v, src_v)
    else
        trafo_ladj = abs(src_v + log(trafo.theta))
        var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
    end
end

const InfiniteDist = Union{
    StandardUvNormal, StandardMvNormal,
    Cauchy, Gumbel, Laplace, Logistic, NoncentralT, Normal, MvNormal, NormalCanon, TDist
}

const PositiveDist = Union{
    BetaPrime, Chi, Chisq, Erlang, Exponential, FDist, Frechet, Gamma, InverseGamma,
    InverseGaussian, Kolmogorov, LogNormal, NoncentralChisq, NoncentralF, Rayleigh, Weibull
}

const UnitDist = Union{
    StandardUvUniform, StandardMvUniform,
    Beta, KSOneSided, NoncentralBeta
}

const UCTruncatedDist = Truncated{<:Distribution{Univariate,Continuous}}

const StdUvDist = Union{StandardUvUniform, StandardUvNormal}
const StdMvDist = Union{StandardMvUniform, StandardMvNormal}

getspace(d::InfiniteDist) = InfiniteSpace()
getspace(d::UnitDist) = UnitSpace()
getspace(d::Distribution) = MixedSpace()

getvartrafo(space::InfiniteSpace, dist::Normal) = IdentityVT(InfiniteSpace(), ScalarShape{Real}())
getvartrafo(space::InfiniteSpace, dist::PositiveDist) = BAT.ScaledLogTrafo(1)
getvartrafo(space::InfiniteSpace, dist::Uniform) = inv(LogisticCDFTrafo(0, 1)) ∘ UniformCDFTrafo(minimum(dist), maximum(dist))
getvartrafo(space::InfiniteSpace, dist::UnitDist) = inv(LogisticCDFTrafo(0, 1))
getvartrafo(space::InfiniteSpace, dist::AbstractMvNormal) = IdentityVT(InfiniteSpace(), ArrayShape{Real}(length(dist)))
# getvartrafo(space::InfiniteSpace, dist::Distributions.AbstractMvLogNormal) = ...