text:
module Franklin

using FranklinTemplates
using FranklinTemplates: filecmp

using Markdown
using Markdown: htmlesc
using Dates
using DelimitedFiles: readdlm
using OrderedCollections
using Pkg
using DocStringExtensions: SIGNATURES, TYPEDEF

import Logging
import LiveServer
import Base.push!
import NodeJS
import Literate
import HTTP
import Random
import ExprTools: splitdef, combinedef
import REPL.REPLCompletions: emoji_symbols

export serve, publish, cleanpull, newsite, optimize, fd2html,
       literate_folder, verify_links, @OUTPUT, get_url

# Extra functions
export lunr

# Convenience functions that get loaded in every module, see eval/module
export fdplotly

# -----------------------------------------------------------------------------
# Legacy with JuDoc
export jd2html # = fd2html

# -----------------------------------------------------------------------------
#
# CONSTANTS
#

const HOME =
    "https://franklinjl.org"
const POINTER_PV =
    "- page variables and {{...}} blocks: $HOME/syntax/page-variables/"
const POINTER_HFUN =
    "- h-functions (hfun) and lx-functions: $HOME/syntax/utils/"
const POINTER_EVAL =
    "- evaluation of Julia code blocks: $HOME/code/"
const POINTER_WORKFLOW =
    "- workflow and folder structure: $HOME/workflow/"

# These are common IP addresses that we can quickly ping to see if the
# user seems online. This is used in `verify_links`. The IPs were
# obtained via `dig www...`; they may change over time; see `check_ping`
# we check in sequence, one should work if the user is online...
const IP_CHECK = (
    "172.217.21.132" => "Google", # google
    "140.82.118.4"   => "GitHub",   # github
    "103.235.46.39"  => "Baidu",  # baidu
    )

"""Big number when we want things to be far."""
const BIG_INT = typemax(Int)

# This dictionary keeps track of a few useful flags that help keept track
# of the context in which particular functions are called. This is not great
# but works fine for what we do here and limits the necessity of passing lots
# of arguments to lots of functions.
const FD_ENV = LittleDict(
    :FULL_PASS     => true,
    :CLEAR         => false,
    :VERB          => false,
    :FINAL_PASS    => false,
    :PRERENDER     => false,
    :NO_FAIL_PRERENDER => true,  # skip prerendering if fails on a page
    :ON_WRITE      => (_, _) -> nothing,
    :FORCE_REEVAL  => false,
    :CUR_PATH      => "",        # complements fd-rpath
    :SOURCE        => "",        # keeps track of the origin of a HTML string
    :OFFSET_LXDEFS => -BIG_INT,  # helps keep track of order in lxcoms/envs
    :DEBUG_MODE    => false,
    :SUPPRESS_ERR  => false,
    :SILENT_MODE   => false,
    :QUIET_TEST    => false,
    :SHOW_WARNINGS => true,     # franklin-specific warnings
    )

# keep track of pages which need to be re-evaluated after the full-pass
# to ensure that their h-fun are working with the fully-defined scope
# (e.g. if need list of all tags)
const DELAYED = Set{String}()

"""Dict to keep track of languages and how comments are indicated and their extensions. This is relevant to allow hiding lines of code. """
const CODE_LANG = LittleDict{String,NTuple{2,String}}(
    # expected most common ones
    "julia"      => (".jl",   "#"),
    "julia-repl" => (".jl",   "#"),
    "python"     => (".py",   "#"),
    "r"          => (".r",    "#"),
    "matlab"     => (".m",    "%"),
    # other ones that could appear (this may get completed over time)
    # note: HTML, Markdown are not here **on purpose** as there can be
    # ambiguities in context.
    "bash"       => (".sh",   "#"),
    "c"          => (".c",    "//"),
    "cpp"        => (".cpp",  "//"),
    "fortran"    => (".f90",  "!"),
    "go"         => (".go",   "//"),
    "javascript" => (".js",   "//"),
    "lua"        => (".lua",  "--"),
    "toml"       => (".toml", "#"),
    "ruby"       => (".ruby", "#"),
    )

# copied from Base/path.jl
if Sys.isunix()
    """Indicator for directory separation on the OS."""
    const PATH_SEP = "/"
elseif Sys.iswindows()
    const PATH_SEP = "\\"
else
    error("Unhandled OS")
end

"""Type of the containers for page variables (local and global)."""
const PageVars = LittleDict{String,Pair}

"""Shorter name for a type that we use everywhere"""
const AS = Union{String,SubString{String}}

"""Convenience constants for an automatic message to add to code files."""
const MESSAGE_FILE_GEN     = "This file was generated, do not modify it."
const MESSAGE_FILE_GEN_LIT = "# $MESSAGE_FILE_GEN\n\n"
const MESSAGE_FILE_GEN_FMD = "# $MESSAGE_FILE_GEN # hide\n"

# -----------------------------------------------------------------------------

"""Debugging mode"""
const LOGGING = Ref(false)

function logger(t)
    LOGGING[] || return nothing
    printstyled("LOG: ", color=:yellow, bold=true)
    printstyled(rpad(t[1], 20), color=:blue)
    println(escape_string(t[2]))
    return nothing
end

# -----------------------------------------------------------------------------

include("build.jl") # check if user has Node/minify

include("regexes.jl")

# UTILS
include("utils/warnings.jl")
include("utils/errors.jl")
include("utils/paths.jl")
include("utils/vars.jl")
include("utils/misc.jl")
include("utils/html.jl")

# PARSING
include("parser/tokens.jl")
include("parser/ocblocks.jl")
# > markdown
include("parser/markdown/tokens.jl")
include("parser/markdown/indent.jl")
include("parser/markdown/validate.jl")
# > latex
include("parser/latex/tokens.jl")
include("parser/latex/blocks.jl")
# > html
include("parser/html/tokens.jl")
include("parser/html/blocks.jl")

# EVAL
include("eval/module.jl")
include("eval/run.jl")
include("eval/codeblock.jl")
include("eval/io.jl")
include("eval/literate.jl")

# CONVERSION
# > markdown
include("converter/markdown/blocks.jl")
include("converter/markdown/utils.jl")
include("converter/markdown/mddefs.jl")
include("converter/markdown/tags.jl")
include("converter/markdown/md.jl")
# > latex
include("converter/latex/latex.jl")
include("converter/latex/objects.jl")
include("converter/latex/hyperrefs.jl")
include("converter/latex/io.jl")
# > html
include("converter/html/functions.jl")
include("converter/html/html.jl")
include("converter/html/blocks.jl")
include("converter/html/link_fixer.jl")
include("converter/html/prerender.jl")

# FILE AND DIR MANAGEMENT
include("manager/rss_generator.jl")
include("manager/sitemap_generator.jl")
include("manager/write_page.jl")
include("manager/dir_utils.jl")
include("manager/file_utils.jl")
include("manager/franklin.jl")
include("manager/extras.jl")
include("manager/post_processing.jl")

end # module


text:
"""
`module  JAC.PlasmaShift`
    ... a submodel of JAC that contains all methods for computing plasma shifts and related properties between for some level or
        between some initial and final-state multiplets.
"""
module PlasmaShift

    using Printf, ..Basics, ..Defaults, ..Nuclear, ..ManyElectron, ..Radial, ..TableStrings

    #============================
    """
    `@enum   PlasmaModel`  ... defines a enumeration for the (allowed) plasma models.

        + NoPlasmaModel      ... No plasma model defined.
        + DebyeHueckel       ... Debye-Hueckel plasma model.
        + IonSphere          ... Ion-sphere (not yet supported).
    """
    @enum   PlasmaModel    NoPlasmaModel    DebyeHueckel    DebeyBox    IonSphere

    """
    `PlasmaModel(sa::String)`  ... constructor for a given String.
    """
    function PlasmaModel(sa::String)
        if       sa == "No Model"              wa = NoPlasmaModel
        elseif   sa == "Debye-Hueckel"         wa = DebyeHueckel
        elseif   sa == "DebeyBox"              wa = DebeyBox
        elseif   sa == "IonSphere"             wa = IonSphere
        else     error("stop a")
        end

        return( UseGauge(wa) )
    end

    # `Base.show(io::IO, model::PlasmaModel)`  ... prepares a proper printout of the variable model::PlasmaModel.
    function Base.show(io::IO, model::PlasmaModel)
        print(io, string(model) )
    end

    # `Base.string(model::PlasmaModel)`  ... provides a proper printout of the variable model::PlasmaModel.
    function Base.string(model::PlasmaModel)
        if       model == NoPlasmaModel     return("No plasma model")
        elseif   model == DebyeHueckel      return("Debye-Hueckel model")
        elseif   model == DebeyBox          return("Debey-box model")
        elseif   model == IonSphere         return("Ion-sphere model")
        else     error("stop a")
        end
    end          ==========================================#

    """
    `abstract type PlasmaShift.AbstractPlasmaModel`
        ... defines an abstract and a number of singleton types for the the (allowed) plasma models.

        + NoPlasmaModel      ... No plasma model defined.
        + DebyeHueckel       ... Debye-Hueckel plasma model.
        + IonSphere          ... Ion-sphere (not yet supported).
    """
    abstract type  AbstractPlasmaModel                      end
    struct         NoPlasmaModel  <:  AbstractPlasmaModel   end
    struct         DebyeHueckel   <:  AbstractPlasmaModel   end
    struct         IonSphere      <:  AbstractPlasmaModel   end

    # `Base.show(io::IO, model::AbstractPlasmaModel)`  ... prepares a proper printout of the variable model::AbstractPlasmaModel.
    function Base.show(io::IO, model::AbstractPlasmaModel)
        print(io, string(model) )
    end

    # `Base.string(model::AbstractPlasmaModel)`  ... provides a proper printout of the variable model::AbstractPlasmaModel.
    function Base.string(model::AbstractPlasmaModel)
        if       model == NoPlasmaModel()     return("No plasma model")
        elseif   model == DebyeHueckel()      return("Debye-Hueckel model")
        elseif   model == DebeyBox()          return("Debey-box model")
        elseif   model == IonSphere()         return("Ion-sphere model")
        else     error("stop a")
        end
    end

    """
    `struct  PlasmaShift.Settings`  ... defines a type for the details and parameters of computing level energies with plasma interactions.

        + plasmaModel      ::AbstractPlasmaModel        ... Specify a particular plasma model, e.g. ion-sphere, debye.
        + lambdaDebye      ::Float64                    ... The lambda parameter of different plasma models.
        + ionSphereR0      ::Float64                    ... The effective radius of the ion-sphere model.
        + NoBoundElectrons ::Int64                      ... Effective number of bound electrons.
    """
    struct Settings
        plasmaModel        ::AbstractPlasmaModel
        lambdaDebye        ::Float64
        ionSphereR0        ::Float64
        NoBoundElectrons   ::Int64
    end

    """
    `PlasmaShift.Settings()`  ... constructor for a standard instance of PlasmaShift.Settings.
    """
    function Settings()
        Settings(DebyeHueckel(), 0.25, 0., 0)
    end

    # `Base.show(io::IO, settings::PlasmaShift.Settings)`  ... prepares a proper printout of the settings::PlasmaShift.Settings.
    function Base.show(io::IO, settings::PlasmaShift.Settings)
        println(io, "plasmaModel:            $(settings.plasmaModel)  ")
        println(io, "lambdaDebye:            $(settings.lambdaDebye)  ")
        println(io, "ionSphereR0:            $(settings.ionSphereR0)  ")
        println(io, "NoBoundElectrons:       $(settings.NoBoundElectrons)  ")
    end

    """
    `struct  PlasmaShift.AugerSettings`
        ... defines a type for the details and parameters of computing Auger rates with plasma interactions.

        + plasmaModel            ::AbstractPlasmaModel           ... Specify a particular plasma model, e.g. ion-sphere, debye.
        + lambdaDebye            ::Float64                       ... The lambda parameter of different plasma models.
        + ionSphereR0            ::Float64                       ... The effective radius of the ion-sphere model.
        + NoBoundElectrons       ::Int64                         ... Effective number of bound electrons.
        + printBeforeComputation ::Bool                          ... True, if all energies and lines are printed before their evaluation.
        + selectLines            ::Bool                          ... True, if lines are selected individually for the computations.
        + selectedLines          ::Array{Tuple{Int64,Int64},1}   ... List of lines, given by tupels (inital-level, final-level).
    """
    struct AugerSettings
        plasmaModel              ::AbstractPlasmaModel
        lambdaDebye              ::Float64
        ionSphereR0              ::Float64
        NoBoundElectrons         ::Int64
        printBeforeComputation   ::Bool
        selectLines              ::Bool
        selectedLines            ::Array{Tuple{Int64,Int64},1}
    end

    """
    `PlasmaShift.AugerSettings()`  ... constructor for a standard instance of PlasmaShift.AugerSettings.
    """
    function AugerSettings()
        AugerSettings(DebyeHueckel(), 0.25, 0., 0, true, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::PlasmaShift.AugerSettings)`  ... prepares a proper printout of the settings::PlasmaShift.AugerSettings.
    function Base.show(io::IO, settings::PlasmaShift.AugerSettings)
        println(io, "plasmaModel:             $(settings.plasmaModel)  ")
        println(io, "lambdaDebye:             $(settings.lambdaDebye)  ")
        println(io, "ionSphereR0:             $(settings.ionSphereR0)  ")
        println(io, "NoBoundElectrons:        $(settings.NoBoundElectrons)  ")
        println(io, "printBeforeComputation:  $(settings.printBeforeComputation)  ")
        println(io, "selectLines:             $(settings.selectLines)  ")
        println(io, "selectedLines:           $(settings.selectedLines)  ")
    end

    """
    `struct  PlasmaShift.PhotoSettings`
        ... defines a type for the details and parameters of computing photoionization rates with plasma interactions.

        + plasmaModel            ::AbstractPlasmaModel           ... Specify a particular plasma model, e.g. ion-sphere, debye.
        + lambdaDebye            ::Float64                       ... The lambda parameter of different plasma models.
        + ionSphereR0            ::Float64                       ... The effective radius of the ion-sphere model.
        + NoBoundElectrons       ::Int64                         ... Effective number of bound electrons.
        + multipoles             ::Array{Basics.EmMultipole}     ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                 ::Array{Basics.UseGauge}        ... Specifies the gauges to be included into the computations.
        + photonEnergies         ::Array{Float64,1}              ... List of photon energies.
        + printBeforeComputation ::Bool                          ... True, if all energies and lines are printed before their evaluation.
        + selectLines            ::Bool                          ... True, if lines are selected individually for the computations.
        + selectedLines          ::Array{Tuple{Int64,Int64},1}   ... List of lines, given by tupels (inital-level, final-level).
    """
    struct PhotoSettings
        plasmaModel              ::AbstractPlasmaModel
        lambdaDebye              ::Float64
        ionSphereR0              ::Float64
        NoBoundElectrons         ::Int64
        multipoles               ::Array{Basics.EmMultipole}
        gauges                   ::Array{Basics.UseGauge}
        photonEnergies           ::Array{Float64,1}
        printBeforeComputation   ::Bool
        selectLines              ::Bool
        selectedLines            ::Array{Tuple{Int64,Int64},1}
    end

    """
    `PlasmaShift.PhotoSettings()`  ... constructor for a standard instance of PlasmaShift.PhotoSettings.
    """
    function PhotoSettings()
        PhotoSettings(DebyeHueckel(), 0.25, 0., 0, [E1], [Basics.UseCoulomb], Float64[], true, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::PlasmaShift.PhotoSettings)`  ... prepares a proper printout of the settings::PlasmaShift.PhotoSettings.
    function Base.show(io::IO, settings::PlasmaShift.PhotoSettings)
        println(io, "plasmaModel:             $(settings.plasmaModel)  ")
        println(io, "lambdaDebye:             $(settings.lambdaDebye)  ")
        println(io, "ionSphereR0:             $(settings.ionSphereR0)  ")
        println(io, "NoBoundElectrons:        $(settings.NoBoundElectrons)  ")
        println(io, "multipoles:              $(settings.multipoles)  ")
        println(io, "gauges:                  $(settings.gauges)  ")
        println(io, "photonEnergies:          $(settings.photonEnergies)  ")
        println(io, "printBeforeComputation:  $(settings.printBeforeComputation)  ")
        println(io, "selectLines:             $(settings.selectLines)  ")
        println(io, "selectedLines:           $(settings.selectedLines)  ")
    end

    """
    `PlasmaShift.computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                     asfSettings::AsfSettings, settings::PlasmaShift.Settings; output=true)`
        ... to compute a new multiplet from the plasma-modified Hamiltonian matrix as specified by the given
            settings. The diagonalization of the (plasma-modified) Hamiltonian matrix follows the asfSettings that
            were applied before for the computation of the multiplet; warnings are issued if features were selected
            that are not supported for plasma calculations, such as the Breit interaction, QED shifts and others.
            The energies and plasma shifts (with regard to the unperturbed multiplet) are printed to screen
            and into the summary file. The generated (plasma) multiplet is returned if output=true and nothing otherwise,
    """
    function computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                             asfSettings::AsfSettings, settings::PlasmaShift.Settings; output=true)
        println("")
        printstyled("PlasmaShift.computeOutcomes(): The computation of plasma-shifts starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------------------------------------- \n", color=:light_green)
        #
        basis        = multiplet.levels[1].basis
        newMultiplet = Basics.perform("computation: CI for plasma", basis, nm, grid, asfSettings, settings)
        # Print all results to screen
        PlasmaShift.displayResults(stdout, multiplet, newMultiplet, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    PlasmaShift.displayResults(iostream, multiplet, newMultiplet, settings)   end
        #
        if    output    return( newMultiplet )
        else            return( nothing )
        end
    end

    """
    `PlasmaShift.displayResults(stream::IO, multiplet::Multiplet, newMultiplet::Multiplet, settings::PlasmaShift.Settings)`
        ... to display the energies, M_ms and F-parameters, etc. for the  selected levels. A neat table is printed but nothing
            is returned otherwise.
    """
    function  displayResults(stream::IO, multiplet::Multiplet, pMultiplet::Multiplet, settings::PlasmaShift.Settings)
        println(stream, " ")
        println(stream, " ")
        println(stream, "  Plasma shifts for $(settings.plasmaModel):")

        # Write out the relevant plasma parameters
        if settings.plasmaModel == PlasmaShift.DebyeHueckel()
            println(stream, "\n     + lambda = $(settings.lambdaDebye)")
            println(stream,   "     + Plasma screening included perturbatively in CI matrix but not in SCF field.")
        else  error("Unsupported plasma model = $(plasmaSettings.plasmaModel)")
        end

        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(64))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(18, "Energy w/o plasma"; na=4)
        sb = sb * TableStrings.center(18, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(14, "Delta E";     na=4)
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(64))
        #
        for  i  in  1:length(multiplet.levels)
            sa  = "  ";    sym = LevelSymmetry( multiplet.levels[i].J, multiplet.levels[i].parity);
                        newsym = LevelSymmetry( pMultiplet.levels[i].J, pMultiplet.levels[i].parity)
            sa = sa * TableStrings.center(10, TableStrings.level(multiplet.levels[i].index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=6)
            energy = multiplet.levels[i].energy
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))              * "     "
            if  sym == newsym   deltaE = pMultiplet.levels[i].energy - energy;
                                sa     = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", deltaE))
            else                sa     = sa * "Level crossing."
            end
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(64), "\n\n")
        #
        return( nothing )
    end

end # module

ext: jl


text:
"""Implementation of propagation via direct matrix exponentiation."""

using LinearAlgebra
import StaticArrays

"""
```julia
    ExpPropWrk(v0)
```

Workspace for propagation via direct matrix exponentiation.

Initializes the workspace for the propagation of a vector `v0`
"""
struct ExpPropWrk{T}
    v :: T
    function ExpPropWrk(v0::T) where T
        new{T}(similar(v0))
    end
end

"""
```julia
expprop!(Ψ, H, dt, wrk; func=(H_dt -> exp(-1im * H_dt)))
```
Evaluate `Ψ = func(H*dt) Ψ` by directly evaluating `U = func(H*dt)`, i.e. by
matrix exponentiation for the default `func`, and then multiplying `U` and
`Ψ` in-place with `mul!`.

The workspace `wrk` must be initialized with [`ExpPropWrk`](@ref) to provide
storage for a temporary state.
"""
function expprop!(Ψ, H, dt, wrk; kwargs...)
    func = get(kwargs, :func, H_dt -> exp(-1im * H_dt))
    copyto!(wrk.v, Ψ)
    U = func(H*dt)
    mul!(Ψ, U, wrk.v)
end

function expprop(Ψ::StaticArrays.SVector, H, dt, wrk; kwargs...)
    func = get(kwargs, :func, H_dt -> exp(-1im * H_dt))
    return func(H*dt) * Ψ
end

function expprop(Ψ, H, dt, wrk; kwargs...)
    func = get(kwargs, :func, H_dt -> exp(-1im * H_dt))
    return func(H*dt) * Ψ
end


text:
using Test
using FakeServer
using HTTP

using Logging
global_logger(ConsoleLogger(stderr,Logging.Debug))
disable_logging(Logging.LogLevel(-10000))
HTTP.DEBUG_LEVEL[]=10

@testset "all tests" begin
    include("fakeserver_test.jl")

end


text:
"""
    rasterzones!(GI::GItype, shapes::Vector{GMTdataset}, fun::Function)

Apply a unidimensional function `fun` to to the elements of the grid or image `GI` that lie inside the polygons
of the GMTdataset `shapes`. The `GI` array is modified in place.

### Parameters
- `GI`: A grid (GMTgrid) or image (GMTimage) type that will be modified by applying `fun` to the elements that
   fall inside the polygons of `shapes`.
- `shapes`: A vector of GMTdataset containing the polygons inside which the elements if `GI` will be assigned
   a single value obtained by applying the function `fun`.
- `fun`: A unidemensional function name used to compute the contant value for the `GI` elements that fall
   inside each of the polygons of `shapes`.

See also: ``colorzones!``

### Returns
It does't return anything but the input `GI` is modified.

## Example
    Take the Peaks grid and replace the elements that fall inside a triangle at the center by their average.

	G = GMT.peaks();
	D = mat2ds([-1 -1; 0 1; 1 -1; -1 -1]);
	rasterzones!(G, D, mean)

"""
function rasterzones!(GI::GItype, shapes::Vector{GMTdataset}, fun::Function)

	function within(bbox_p, bbox_R)		# Check if the polygon BB is contained inside the image's region.
		bbox_p[1] >= bbox_R[1] && bbox_p[2] <= bbox_R[2] && bbox_p[3] >= bbox_R[3] && bbox_p[4] <= bbox_R[4]
	end
	function maskit(GI_, mask, band)
		if (band > 0)
			r = (eltype(GI_) <: Integer) ? round(eltype(GI_), fun(GI_[mask, band])) : fun(GI_[mask, band])
		else
			r = (eltype(GI_) <: Integer) ? round(eltype(GI_), fun(GI_[mask])) : fun(GI_[mask])
		end
	end

	# GDAL always returns TRB, so if GI has a different one, must convert. Also assume that an empty layout <=> BCB
	layout = startswith(GI.layout, "TR") ? "" : (GI.layout == "" ? "BCB" : GI.layout)
	row_dim, col_dim = (GI.layout == "" || GI.layout[2] == 'C') ? (1,2) : (2,1)	# If RowMajor the array is transposed

	for k = 1:length(shapes)
		!within(shapes[k].bbox, GI.range) && continue		# Catch any exterior polygon before it errors
		_GI, pix_x, pix_y = GMT.crop(GI, region=shapes[k])

		mask = reinterpret(Bool, gdalrasterize(shapes[k],
			["-of", "MEM", "-ts","$(size(_GI, col_dim))","$(size(_GI, row_dim))", "-burn", "1", "-ot", "Byte"], layout=layout))
		(!any(mask)) && continue		# If mask is all falses stop before it errors

		if (ndims(_GI) == 2)
			_GI[mask] .= maskit(_GI, mask, 0)
			GI[pix_y[1]:pix_y[2], pix_x[1]:pix_x[2]] = _GI
		else
			for n = 1:3
				_GI[mask, n] .= maskit(_GI, mask, n)
				GI[pix_y[1]:pix_y[2], pix_x[1]:pix_x[2], n] = _GI[:,:,n]
			end
		end
	end
	nothing
end

# ---------------------------------------------------------------------------------------------------
"""
    colorzones!(shapes::Vector{GMTdataset}[, fun::Function]; img::GMTimage=nothing, url::AbstractString="", layer=0, pixelsize::Int=0, append::Bool=true)

Paint the polygons in the `shapes` with the average color that those polygons ocupy in the `img` image.
When the `shapes` are plotted the resulting image looks like a choropleth map. Alternatively, instead of
transmitting the `img` image one can provide a WMS URL, `layer` number and `pixelsize` to download images from
the Web Map Server service, covering each one the bounding box of each of the `shapes` polygons. This option
is much slower than passing the entire image at once but consumes much less memory. Important when the total area
is large (think Russia size) because even at a moderately resultion it can imply downloading a huge file.

### Parameters
- `shapes`: A vector of GMTdataset containing the polygons to be painted. This container will be changed
   in the sense that the ``header`` field of each polygon (a GMTdataset) will be appended with the fill
   color (but see also the `append` option that controls how this change takes place.)
- `fun`: By default the average color is obtained by taking the square root of the average of squares of each
   of the three (RGB) bands (or just one for grayscale images). Give the name of another function to replace
   this default. For example ``median`` will assign the color by computing the median of each component inside
   the polygon area in question.
- `img`: the image from which the stats (`fun`) of each color for each polygon will be computed.
- `url`: In case the `img` option is not used, pass the Web Map Server URL (see the ``wmsinfo`` and ``wmsread`` functions)
   from where the the images covering the BoundingBox of each polygon will be downloaded. Warning, this is a much slower
   method but potentially usefull when the images to download risk to be very big.
- `layer`: When the `url` option is used this one becomes mandatory and represents the layer number or layer name of
   interest from those provided by the WMS service. That is, both of these forms are allowed: `layer=3` or
   `layer="Invented layer name"`
- `pixelsize`: Sets the requested cell size in meters [default]. Use a string appended with a 'd'
   (e.g. `resolution="0.001d"`) if the resolution is given in degrees. This way works only when the layer is in geogs.
- `append`: By default the color assignement to each of the polygons in the `shapes` vector is achieved by
   appending the fill color to the possibly existing header field. Running the ``colorzones!`` command more than once
   keeps adding (now ignored, because only the first is used) colors. Setting `append=false` forces rewriting
   the header field at each run and hence the assigned color is always the one used (but the previous header is cleared out).

See also: ``rasterzones!``

### Returns
It does't return anything but the input `shapes` is modified.

## Example
    Read the 2020 Sentinel2 Earth color for Portugal at 100 m resolution. Load the administrative
	regions and compute their median colors.

    wms = wmsinfo("http://tiles.maps.eox.at/wms?");
    img = wmsread(wms, layer=3, region=(-9.6,-6,36.9,42.2), pixelsize=100);
	Pt = gmtread("C:\\programs\\compa_libs\\covid19pt\\extra\\mapas\\concelhos\\concelhos.shp");
	colorzones!(Pt, median, img=img);
	imshow(Pt, proj=:guess)
"""
colorzones!(shapes::Vector{GMTdataset}; img::GMTimage=nothing, url::AbstractString="", layer=0, pixelsize::Int=0, append::Bool=true) =
	colorzones!(shapes, meansqrt; img=img, url=url, layer=layer, pixelsize=pixelsize, append=append)
function colorzones!(shapes::Vector{GMTdataset}, fun::Function; img::GMTimage=nothing, url::AbstractString="", layer=0, pixelsize::Int=0, append::Bool=true)

	(img === nothing && url == "") && error("Must either pass a grid/image or a WMS URL.")
	(img !== nothing && url != "") && error("Make up your mind. Pass ONLY one of grid/image or a WMS URL.")
	(url != "" && (layer == 0 || pixelsize == 0)) && error("When passing a WMS URL must also provide the layer and the pixelsize.")

	function within(k)		# Check if the polygon BB is contained inside the image's region.
		shapes[k].bbox[1] >= img.range[1] && shapes[k].bbox[2] <= img.range[2] && shapes[k].bbox[3] >= img.range[3] && shapes[k].bbox[4] <= img.range[4]
	end

	(url != "") && (wms = wmsinfo(url))
	(url != "") && (layer_n = get_layer_number(wms, layer))
	layout = ""		# GDAL always returns TRB, so if img has a different one, we must convert (arg in gdalrasterize).
	(img !== nothing && !startswith(img.layout, "TR")) && (layout = img.layout)
	row_dim, col_dim = (img.layout == "" || img.layout[2] == 'C') ? (1,2) : (2,1)	# If RowMajor the array is transposed

	for k = 1:length(shapes)
		!within(k) && continue				# Catch any exterior polygon before it errors
		_img = (url != "") ? wmsread(wms, layer=layer_n, region=shapes[k], pixelsize=pixelsize) : GMT.crop(img, region=shapes[k])[1]

		mk = gdalrasterize(shapes[k],
			["-of", "MEM", "-ts","$(size(_img, col_dim))","$(size(_img, row_dim))", "-burn", "1", "-ot", "Byte"], layout=layout)
		mask = reinterpret(Bool, mk)
		(!any(mask)) && continue		# If mask is all falses stop before it errors

		opt_G = (ndims(_img) >= 3) ? @sprintf(" -G%.0f/%.0f/%.0f", fun(_img[mask,1]), fun(_img[mask,2]), fun(_img[mask,3])) :
		         (c = fun(_img[mask,1]); @sprintf(" -G%.0f/%.0f/%.0f", c,c,c))
		shapes[k].header = (append) ? shapes[k].header * opt_G : opt_G
		(url != "" && rem(k, 10) == 0) && (println("Done ",k, " of ", length(shapes)); print(stdout, "\e[", 1, "A", "\e[1G"))
	end
	shapes
end

function meansqrt(x)
	sum = 0.0
	@inbounds for k = 1:length(x)
		sum += Float64(x[k]) * x[k]
	end
	sqrt(sum/length(x))
end


text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "output tests for latin" begin
    data_path = joinpath(@__DIR__, "data", "latin_mini.csv")
    output_dir_path = joinpath(@__DIR__, "out_data", "latin_mini.csv")

    JudiLing.lpo_cv_split(10, data_path)
    JudiLing.loo_cv_split(data_path)

    JudiLing.train_val_random_split(
        data_path,
        output_dir_path,
        "latin",
        train_sample_size = 100,
        val_sample_size = 10,
        )

    JudiLing.train_val_random_split(
        data_path,
        output_dir_path,
        "latin",
        val_ratio = 0.1,
        )

    JudiLing.train_val_careful_split(
        data_path,
        output_dir_path,
        "latin",
        ["Lexeme", "Person", "Number", "Tense", "Voice", "Mood"],
        train_sample_size = 100,
        val_sample_size = 10,
        n_grams_target_col = :Word,
        n_grams_tokenized = false,
        n_grams_sep_token = nothing,
        grams = 3,
        n_grams_keep_sep = false,
        start_end_token = "#",
        )

    JudiLing.train_val_careful_split(
        data_path,
        output_dir_path,
        "latin",
        ["Lexeme", "Person", "Number", "Tense", "Voice", "Mood"],
        val_ratio = 0.1,
        n_grams_target_col = :Word,
        n_grams_tokenized = false,
        n_grams_sep_token = nothing,
        grams = 3,
        n_grams_keep_sep = false,
        start_end_token = "#",
        )

    rm(output_dir_path, force = true, recursive = true)
end


text:
module SocialSim

using Dates
using Random
using LightGraphs
using DataFrames
using CSV
using Dates
using GraphIO
using Plots
using StatsBase
using ProgressMeter

export run_sim
export sweep_temp

const J = 1

include("init.jl")
include("analysis.jl")
include("dynamics.jl")
include("storage.jl")

#Runs the simulation at a determined temperature for a determined number of steps
function run_sim(T,population,steps; exports_number=10)

    #Random.seed!(1234)

    #Get time of simulation and prepare folders
    exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")
    mkdir("Data/$exportTime")
    mkdir("Data/$exportTime/Network")

    #Initializing network and data frame
    data = createDataFrame()
    network, nodes = initNetwork(data, population)

    #Export network at each nx steps
    nx = div(steps, exports_number)

    #Executing steps
    @showprogress 5 "Computing..." for i in 1:steps
        ni = rand(1:population)
        procedure2(nodes, network, data, ID=ni, N=population, T=T)

        if mod(i, nx) == 0
            exportNetwork(exportTime, network, nodes, i)
        end
    end

    #Export data to a new folder
    exportData(exportTime, data, network, nodes)
    plotAnalysis(steps, exportTime, data, nodes, network)
end

function sweep_temp(T_f,T_step,population)

    #Random.seed!(1234)

    #Get time of simulation and prepare folders
    exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")
    mkdir("Data/$exportTime")

    #Initializing network and data frame
    data = createDataFrame()
    network, nodes = initNetwork(data, population)

    #Sweep temp
    equilibriumTime = []
    @showprogress 5 "Computing..." for T in 1:T_step:T_f
        counter = 0
        mkdir("Data/$exportTime/T=$T")
        mkdir("Data/$exportTime/T=$T/Network")

        while checkEquilibrium(data) == false
            ni = rand(1:population)
            procedure2(nodes, network, data, ID=ni, N=population, T=T)
            counter += 1
        end

        #Export data to a new folder
        exportData("$exportTime/T=$T", data, network, nodes)
        exportNetwork("$exportTime/T=$T", network, nodes, i)
        plotAnalysis(steps, "$exportTime/T=$T", data, nodes, network)

        push!(equilibriumTime, counter)
        println("Simulation for T=$T done")
    end

    #Export equilibriumTime
    exportEquilibriumLog(T_step, equilibriumTime, exportTime)
    plotEquilibriumTime=plot(1:T_step:T_f, equilibriumTime, legend=false, size(1200,800))
    xlabel!("Temperature")
    ylabel!("Equilibrium Time")
    png("Data/$dir/Equilibrium_time")
end

#Runs the simulation at T for a determined time interval (in hours)
#=function set_run(T,population,duration; exports_number=10)

    #Get time of simulation and prepare folders
    exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")
    mkdir("Data/$exportTime")
    mkdir("Data/$exportTime/Network")

    #Initializing network and data frame
    data = createDataFrame()
    network, nodes = initNetwork(data, population)

    #Export network at each nx steps
    nx = div(steps, exports_number)

    timelimit = Dates.now() + Dates.Minute(duration)

    #Executing steps
    while Dates.now() < timelimit #&& length(Data.E) <= 2000000
        ni = rand(1:population)
        procedure2(nodes, network, data, ID=ni, N=population, T=T)

        if mod(i, nx) == 0
            exportNetwork(exportTime, network, nodes, i)
        end
    end

    #Export data to a new folder
    exportData(exportTime)
    plotAnalysis(length(data.E)-1, exportTime)
end=#

end

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
# we don't introduce a new module here
# for convenience
using Configurations
using .AST

include("types.jl")
include("cast.jl")
include("markdown.jl")
include("utils.jl")


text:
using Print

# strings are defined with double quotes, can contain unicode
s1 = "Greek letter α,β,γ"
println(s1)

# chars are defined with single quotes
c1 = 'a'
println(c1)

# ascii value
println("Int('α') == ", Int('α'))

# upper case or lower case:
s1_caps = uppercase(s1)
s1_lower = lowercase(s1)

# substring
println(s1[1:8])
#> Greek le

# interpolation
a = "Have a"
b = "day"
println("$a nice $b")

# extend
println("1 * 2 = $(1 * 2)")

# concatenated with "*" and string()
println("Hello " * "World!")
println(string("Hello ", "World!"))

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia


text:
#=
block_boundary_experiment_profile:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-04-26
=#

# include("../src/LarSurf.jl")
tim = time()
using Revise
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pandas
using Seaborn
using Dates
using Logging

using Profile
using ProfileView

# using Logger
# logger = SimpleLogger(stdout, Logging.Debug)
logger = SimpleLogger(stdout, Logging.Info)
old_logger = global_logger(logger)

# Logging.min_enabled_level

data_file = "data.csv"
fn_profile_txt = "exp_surface_per_block_profile.txt"

tim_prev = tim
tim = time()
println("using done in: ", tim - tim_prev)

## Read data from file
# xystep = 1
# zstep = 1
# threshold = 4000;
# pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")

threshold = 10
# xystep = 50
# zstep = 30

# xystep = 25
# zstep = 15

# xystep = 20
# zstep = 10

# xystep = 10
# zstep = 5

# xystep = 4
# zstep = 2

xystep = 8
zstep = 4
pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.1/MASKS_DICOM/liver")

datap = Pio3d.read3d(pth);
#
data3d_full = datap["data3d"]
println("orig size: ", size(data3d_full))
data3d = data3d_full[1:zstep:end, 1:xystep:end, 1:xystep:end];

data_size = LarSurf.size_as_array(size(data3d))
println("data size: ", data_size)
segmentation = data3d .> threshold;

tim_prev = tim
tim = time()
println("data read complete in time: ", tim - tim_prev)

# Run once to force compilation
## Generate data
# segmentation = LarSurf.generate_slope([9,10,11])
block_size = [8,8,8]
function time_profile()
    larmodel = LarSurf.get_surface_grid_per_block(segmentation, block_size)
    # bigVV, bigEV, bigFV, bigCV = tmodel
end

#compile
LarSurf.set_param(boundary_allow_read_files=false)
LarSurf.set_param(boundary_allow_write_files=false)
LarSurf.set_param(boundary_allow_memory=false)
time_profile()
println("==========")
# LarSurf._boundary3_storage = Dict()
LarSurf.set_param(boundary_allow_read_files=false)
LarSurf.set_param(boundary_allow_write_files=false)
LarSurf.set_param(boundary_allow_memory=true)
# delete memory
LarSurf.reset(boundary_storage=true)

Profile.init(n = 10^9, delay=0.01)
@profile time_profile()
open(fn_profile_txt, "w") do s
    Profile.print(IOContext(s, :displaysize => (24, 500)))
end
Profile.print(format=:flat)
ProfileView.view()

# ProfileView.svgwrite("profile_results2.svg")


text:
int_rules_4_1_7 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.7*(d*trig)^m*(a+b*(c*sin)^n)^p =#
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (((4 * ~A * (2 * ~a + ~b) + ~B * (4 * ~a + 3 * ~b)) * ~x) / 8 - ((4 * ~A * ~b + ~B * (4 * ~a + 3 * ~b)) * cos(~e + ~f * ~x) * sin(~e + ~f * ~x)) / (8 * ~f)) - (~b * ~B * cos(~e + ~f * ~x) * sin(~e + ~f * ~x) ^ 3) / (4 * ~f) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~B) * cos(~e + ~f * ~x) * sin(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ ~p) / (2 * ~f * (~p + 1)) + (1 / (2 * (~p + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p - 1) * Simp(~a * ~B + 2 * ~a * ~A * (~p + 1) + (2 * ~A * ~b * (~p + 1) + ~B * (~b + 2 * ~a * ~p + 2 * ~b * ~p)) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x) && GtQ(~p, 0)
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x) ^ 2) / (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~B * ~x) / ~b + ((~A * ~b - ~a * ~B) / ~b) * Antiderivative(1 / (~a + ~b * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x)
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x) ^ 2) / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~B / ~b) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x) ^ 2), ~x) + ((~A * ~b - ~a * ~B) / ~b) * Antiderivative(1 / sqrt(~a + ~b * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-((~A * ~b - ~a * ~B)) * cos(~e + ~f * ~x) * sin(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p + 1)) / (2 * ~a * ~f * (~a + ~b) * (~p + 1)) - (1 / (2 * ~a * (~a + ~b) * (~p + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p + 1) * Simp((~a * ~B - ~A * (2 * ~a * (~p + 1) + ~b * (2 * ~p + 3))) + 2 * (~A * ~b - ~a * ~B) * (~p + 2) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x) && (LtQ(~p, -1) && NeQ(~a + ~b, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((ff * (~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ ~p * (sec(~e + ~f * ~x) ^ 2) ^ ~p) / (~f * (~a + (~a + ~b) * tan(~e + ~f * ~x) ^ 2) ^ ~p)) * Subst(Antiderivative(((~a + (~a + ~b) * ff ^ 2 * (~x) ^ 2) ^ ~p * (~A + (~A + ~B) * ff ^ 2 * (~x) ^ 2)) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~p + 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x) && Not(IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => (~a) ^ ~p * Antiderivative(ActivateTrig(~u * cos(~e + ~f * ~x) ^ (2 * ~p)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (EqQ(~a + ~b, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => Antiderivative(ActivateTrig(~u * (~a * cos(~e + ~f * ~x) ^ 2) ^ ~p), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && EqQ(~a + ~b, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (sqrt(~a) / ~f) * Elliptic.E(~e + ~f * ~x, -(~b) / ~a) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && GtQ(~a, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (sqrt(~a + ~b * sin(~e + ~f * ~x) ^ 2) / sqrt(1 + (~b * sin(~e + ~f * ~x) ^ 2) / ~a)) * Antiderivative(sqrt(1 + (~b * sin(~e + ~f * ~x) ^ 2) / ~a), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && Not(GtQ(~a, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ 2, ~x) => (((8 * (~a) ^ 2 + 8 * ~a * ~b + 3 * (~b) ^ 2) * ~x) / 8 - (~b * (8 * ~a + 3 * ~b) * cos(~e + ~f * ~x) * sin(~e + ~f * ~x)) / (8 * ~f)) - ((~b) ^ 2 * cos(~e + ~f * ~x) * sin(~e + ~f * ~x) ^ 3) / (4 * ~f) <-- FreeQ([~a, ~b, ~e, ~f], ~x)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => (-(~b) * cos(~e + ~f * ~x) * sin(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p - 1)) / (2 * ~f * ~p) + (1 / (2 * ~p)) * Antiderivative((~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p - 2) * Simp(~a * (~b + 2 * ~a * ~p) + ~b * (2 * ~a + ~b) * (2 * ~p - 1) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (NeQ(~a + ~b, 0) && GtQ(~p, 1))
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(1 / (~a + (~a + ~b) * ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / (sqrt(~a) * ~f)) * Elliptic.F(~e + ~f * ~x, -(~b) / ~a) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && GtQ(~a, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (sqrt(1 + (~b * sin(~e + ~f * ~x) ^ 2) / ~a) / sqrt(~a + ~b * sin(~e + ~f * ~x) ^ 2)) * Antiderivative(1 / sqrt(1 + (~b * sin(~e + ~f * ~x) ^ 2) / ~a), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && Not(GtQ(~a, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => (-(~b) * cos(~e + ~f * ~x) * sin(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p + 1)) / (2 * ~a * ~f * (~p + 1) * (~a + ~b)) + (1 / (2 * ~a * (~p + 1) * (~a + ~b))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x) ^ 2) ^ (~p + 1) * Simp((2 * ~a * (~p + 1) + ~b * (2 * ~p + 3)) - 2 * ~b * (~p + 2) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (NeQ(~a + ~b, 0) && LtQ(~p, -1))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff * sqrt(cos(~e + ~f * ~x) ^ 2)) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative((~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p / sqrt(1 - ff ^ 2 * (~x) ^ 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && Not(IntegerQ(~p))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (-ff / ~f) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * ((~a + ~b) - ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ((~m - 1) / 2)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~a + (~a + ~b) * ff ^ 2 * (~x) ^ 2) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && IntegerQ(~p))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff ^ (~m + 1) * sqrt(cos(~e + ~f * ~x) ^ 2)) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p) / sqrt(1 - ff ^ 2 * (~x) ^ 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~m / 2) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], ((-ff * (~d) ^ (2 * IntPart((~m - 1) / 2) + 1) * (~d * sin(~e + ~f * ~x)) ^ (2 * FracPart((~m - 1) / 2))) / (~f * (sin(~e + ~f * ~x) ^ 2) ^ FracPart((~m - 1) / 2))) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * ((~a + ~b) - ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ((~m - 1) / 2)
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + (~a + ~b) * ff ^ 2 * (~x) ^ 2) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && IntegerQ(~p))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff * sqrt(cos(~e + ~f * ~x) ^ 2)) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~m / 2) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~(d') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff * (~d) ^ (2 * IntPart((~m - 1) / 2) + 1) * (~d * cos(~e + ~f * ~x)) ^ (2 * FracPart((~m - 1) / 2))) / (~f * (cos(~e + ~f * ~x) ^ 2) ^ FracPart((~m - 1) / 2))) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x) ^ 2, ~x)], (ff ^ ((~m + 1) / 2) / (2 * ~f)) * Subst(Antiderivative(((~x) ^ ((~m - 1) / 2) * (~a + ~b * ff * ~x) ^ ~p) / (1 - ff * ~x) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) ^ 2 / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ((~m - 1) / 2)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(((~d * ff * ~x) ^ ~m * (~a + (~a + ~b) * ff ^ 2 * (~x) ^ 2) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && IntegerQ(~p)
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff ^ (~m + 1) * sqrt(cos(~e + ~f * ~x) ^ 2)) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p) / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~m / 2) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff * (~d * tan(~e + ~f * ~x)) ^ (~m + 1) * (cos(~e + ~f * ~x) ^ 2) ^ ((~m + 1) / 2)) / (~d * ~f * sin(~e + ~f * ~x) ^ (~m + 1))) * Subst(Antiderivative(((ff * ~x) ^ ~m * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p) / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(m') * (~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~d * ff * ~x) ^ ~n * (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~n, ~p], ~x) && IntegerQ((~m - 1) / 2)
    @apply_utils Antiderivative((~(c') * sin(~(e') + ~(f') * ~x)) ^ ~m * sin(~(e') + ~(f') * ~x) ^ ~(n') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (-ff / ~f) * Subst(Antiderivative((~c * ff * ~x) ^ ~m * (1 - ff ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2) * ((~a + ~b) - ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~m, ~p], ~x) && IntegerQ((~n - 1) / 2)
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m * sin(~(e') + ~(f') * ~x) ^ ~n * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff ^ (~n + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~n * (~a + (~a + ~b) * ff ^ 2 * (~x) ^ 2) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ ((~m + ~n) / 2 + ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~n / 2) && IntegerQ(~p)))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m * (~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff * sqrt(cos(~e + ~f * ~x) ^ 2)) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative((~d * ff * ~x) ^ ~n * (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~n, ~p], ~x) && IntegerQ(~m / 2)
    @apply_utils Antiderivative((~(c') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], ((ff * (~c) ^ (2 * IntPart((~m - 1) / 2) + 1) * (~c * cos(~e + ~f * ~x)) ^ (2 * FracPart((~m - 1) / 2))) / (~f * (cos(~e + ~f * ~x) ^ 2) ^ FracPart((~m - 1) / 2))) * Subst(Antiderivative((~d * ff * ~x) ^ ~n * (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => (-(cot(~e + ~f * ~x)) * (~b * sin(~e + ~f * ~x) ^ 2) ^ ~p) / (2 * ~f * ~p) + ((~b * (2 * ~p - 1)) / (2 * ~p)) * Antiderivative((~b * sin(~e + ~f * ~x) ^ 2) ^ (~p - 1), ~x) <-- FreeQ([~b, ~e, ~f], ~x) && (Not(IntegerQ(~p)) && GtQ(~p, 1))
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => (cot(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x) ^ 2) ^ (~p + 1)) / (~b * ~f * (2 * ~p + 1)) + ((2 * (~p + 1)) / (~b * (2 * ~p + 1))) * Antiderivative((~b * sin(~e + ~f * ~x) ^ 2) ^ (~p + 1), ~x) <-- FreeQ([~b, ~e, ~f], ~x) && (Not(IntegerQ(~p)) && LtQ(~p, -1))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x) ^ 2, ~x)], (ff ^ ((~m + 1) / 2) / (2 * ~f)) * Subst(Antiderivative(((~x) ^ ((~m - 1) / 2) * (~b * ff ^ (~n / 2) * (~x) ^ (~n / 2)) ^ ~p) / (1 - ff * ~x) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) ^ 2 / ff)) <-- FreeQ([~b, ~e, ~f, ~p], ~x) && (IntegerQ((~m - 1) / 2) && IntegerQ(~n / 2))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~b * (~c * ff * ~x) ^ ~n) ^ ~p) / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~b, ~c, ~e, ~f, ~n, ~p], ~x) && ILtQ((~m - 1) / 2, 0)
    @apply_utils Antiderivative(~(u') * (~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (((~b * ff ^ ~n) ^ IntPart(~p) * (~b * sin(~e + ~f * ~x) ^ ~n) ^ FracPart(~p)) / (sin(~e + ~f * ~x) / ff) ^ (~n * FracPart(~p))) * Antiderivative(ActivateTrig(~u) * (sin(~e + ~f * ~x) / ff) ^ (~n * ~p), ~x)) <-- FreeQ([~b, ~e, ~f, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (IntegerQ(~n) && (EqQ(~u, 1) || MatchQ(~u, (~(d') * (~trig)(~e + ~f * ~x)) ^ ~(m') <-- FreeQ([d, m], ~x) && MemberQ([sin, cos, tan, cot, sec, csc], trig)))))
    @apply_utils Antiderivative(~(u') * (~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (((~b) ^ IntPart(~p) * (~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ FracPart(~p)) / (~c * sin(~e + ~f * ~x)) ^ (~n * FracPart(~p))) * Antiderivative(ActivateTrig(~u) * (~c * sin(~e + ~f * ~x)) ^ (~n * ~p), ~x) <-- FreeQ([~b, ~c, ~e, ~f, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~n)) && (EqQ(~u, 1) || MatchQ(~u, (~(d') * (~trig)(~e + ~f * ~x)) ^ ~(m') <-- FreeQ([d, m], ~x) && MemberQ([sin, cos, tan, cot, sec, csc], trig)))))
    #= Antiderivative(((~a)+(~b')*sin((~e')+(~f')*(~x))^4)^(~p'),~x) :=  With([ff=FreeFactors(tan(e+f*x),x)], -ff/f*Subst(Antiderivative((a+b+2*a*ff^2*x^2+a*ff^4*x^4)^p/(1+ff^2*x^2)^(2*p+1) ,x),x,cot(e+f*x)/ff)) <-- FreeQ([a,b,e,f],x) && IntegerQ(p) =#
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + 2 * ~a * ff ^ 2 * (~x) ^ 2 + (~a + ~b) * ff ^ 4 * (~x) ^ 4) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2) ^ (2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IntegerQ(~p)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((ff * (~a + ~b * sin(~e + ~f * ~x) ^ 4) ^ ~p * (sec(~e + ~f * ~x) ^ 2) ^ (2 * ~p)) / (~f * (~a + 2 * ~a * tan(~e + ~f * ~x) ^ 2 + (~a + ~b) * tan(~e + ~f * ~x) ^ 4) ^ ~p)) * Subst(Antiderivative((~a + 2 * ~a * ff ^ 2 * (~x) ^ 2 + (~a + ~b) * ff ^ 4 * (~x) ^ 4) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2) ^ (2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ(~p - 1 / 2)
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n), ~x) => Module([k], Dist(2 / (~a * ~n), Sum(Antiderivative(1 / (1 - sin(~e + ~f * ~x) ^ 2 / ((-1) ^ ((4k) / ~n) * Rt(-(~a) / ~b, ~n / 2))), ~x), [k, 1, ~n / 2]), ~x)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IntegerQ(~n / 2)
    #= Antiderivative(((~a)+(~b')*sin((~e')+(~f')*(~x))^(~n))^(~p),~x) :=  With([ff=FreeFactors(tan(e+f*x),x)], -ff/f*Subst(Antiderivative((b+a*(1+ff^2*x^2)^(n/2))^p/(1+ff^2*x^2)^(n*p/2+1),x) ,x,cot(e+f*x)/ff)) <-- FreeQ([a,b,e,f],x) && IntegerQ(n/2) &&  IGtQ(p,0) =#
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~b * ff ^ ~n * (~x) ^ ~n + ~a * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2) ^ ((~n * ~p) / 2 + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~n / 2) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => Antiderivative(ExpandTrig((~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n], ~x) && (IGtQ(~p, 0) || EqQ(~p, -1) && IntegerQ(~n))
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => Unintegrable((~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (-ff / ~f) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (((~a + ~b) - 2 * ~b * ff ^ 2 * (~x) ^ 2) + ~b * ff ^ 4 * (~x) ^ 4) ^ ~p, ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ((~m - 1) / 2)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (-ff / ~f) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * (1 - ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p, ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ((~m - 1) / 2) && IntegerQ(~n / 2))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~a + 2 * ~a * ff ^ 2 * (~x) ^ 2 + (~a + ~b) * ff ^ 4 * (~x) ^ 4) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && IntegerQ(~p))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~a * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~b * ff ^ ~n * (~x) ^ ~n) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~p) / 2 + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~n / 2) && IntegerQ(~p)))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((ff ^ (~m + 1) * (~a + ~b * sin(~e + ~f * ~x) ^ 4) ^ ~p * (sec(~e + ~f * ~x) ^ 2) ^ (2 * ~p)) / (~f * Apart(~a * (1 + tan(~e + ~f * ~x) ^ 2) ^ 2 + ~b * tan(~e + ~f * ~x) ^ 4) ^ ~p)) * Subst(Antiderivative(((~x) ^ ~m * ExpandToSum(~a * (1 + ff ^ 2 * (~x) ^ 2) ^ 2 + ~b * ff ^ 4 * (~x) ^ 4, ~x) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~m / 2) && IntegerQ(~p - 1 / 2))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig(sin(~e + ~f * ~x) ^ ~m * (~a + ~b * sin(~e + ~f * ~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegersQ(~m, ~p) && (EqQ(~n, 4) || (GtQ(~p, 0) || EqQ(~p, -1) && IntegerQ(~n))))
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig((~d * sin(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * sin(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * (~c * ff * ~x) ^ ~n) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && (IntegerQ((~m - 1) / 2) && (EqQ(~n, 4) || (GtQ(~m, 0) || (IGtQ(~p, 0) || IntegersQ(~m, ~p)))))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + 2 * ~a * ff ^ 2 * (~x) ^ 2 + (~a + ~b) * ff ^ 4 * (~x) ^ 4) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && IntegerQ(~p))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~b * ff ^ ~n * (~x) ^ ~n + ~a * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + (~n * ~p) / 2 + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ(~m / 2) && (IntegerQ(~n / 2) && IntegerQ(~p)))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~m / (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n), ~x) => Antiderivative(expand((1 - sin(~e + ~f * ~x) ^ 2) ^ (~m / 2) / (~a + ~b * sin(~e + ~f * ~x) ^ ~n), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IGtQ(~m / 2, 0) && IntegerQ((~n - 1) / 2))
    #= Antiderivative(cos((~e')+(~f')*(~x))^(~m)*((~a)+(~b')*sin((~e')+(~f')*(~x))^(~n))^(~p),~x) :=   Antiderivative(ExpandTrig((1-sin(e+f*x)^2)^(m/2)*(a+b*sin(e+f*x)^n)^p,x),x) <--  FreeQ([a,b,e,f],x) && IntegerQ(m/2) && IntegerQ((n-1)/2) &&  ILtQ(p,-1) && LtQ(m,0) =#
    @apply_utils Antiderivative((~(d') * cos(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig((~d * cos(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(d') * cos(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * cos(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x) ^ 2, ~x)], (ff ^ ((~m + 1) / 2) / (2 * ~f)) * Subst(Antiderivative(((~x) ^ ((~m - 1) / 2) * (~a + ~b * ff ^ (~n / 2) * (~x) ^ (~n / 2)) ^ ~p) / (1 - ff * ~x) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) ^ 2 / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ((~m - 1) / 2) && IntegerQ(~n / 2))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * (~c * ff * ~x) ^ ~n) ^ ~p) / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && ILtQ((~m - 1) / 2, 0)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(((~d * ff * ~x) ^ ~m * ExpandToSum(~a * (1 + ff ^ 2 * (~x) ^ 2) ^ 2 + ~b * ff ^ 4 * (~x) ^ 4, ~x) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && IntegerQ(~p)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ 4) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((ff * (~a + ~b * sin(~e + ~f * ~x) ^ 4) ^ ~p * (sec(~e + ~f * ~x) ^ 2) ^ (2 * ~p)) / (~f * Apart(~a * (1 + tan(~e + ~f * ~x) ^ 2) ^ 2 + ~b * tan(~e + ~f * ~x) ^ 4) ^ ~p)) * Subst(Antiderivative(((~d * ff * ~x) ^ ~m * ExpandToSum(~a * (1 + ff ^ 2 * (~x) ^ 2) ^ 2 + ~b * ff ^ 4 * (~x) ^ 4, ~x) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (2 * ~p + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && IntegerQ(~p - 1 / 2)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~d * ~x) ^ ~m * (~b * ff ^ ~n * (~x) ^ ~n + ~a * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ ((~n * ~p) / 2 + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && (IntegerQ(~n / 2) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig((~d * tan(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * tan(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * cot(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * cot(~e + ~f * ~x)) ^ FracPart(~m) * (tan(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((tan(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(d') * sec(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * sec(~e + ~f * ~x)) ^ FracPart(~m) * (cos(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((cos(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sin(~(e') + ~(f') * ~x) ^ ~(n')) ^ ~(p'), ~x) => (~d) ^ (~n * ~p) * Antiderivative((~d * csc(~e + ~f * ~x)) ^ (~m - ~n * ~p) * (~b + ~a * csc(~e + ~f * ~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~m)) && IntegersQ(~n, ~p))
    @apply_utils Antiderivative((~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * csc(~e + ~f * ~x)) ^ FracPart(~m) * (sin(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((sin(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * sin(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x) + ~(d') * cos(~(e') + ~(f') * ~x)) ^ 2) ^ ~p, ~x) => Antiderivative((~a + ~b * (sqrt((~c) ^ 2 + (~d) ^ 2) * sin(atan(~c, ~d) + ~e + ~f * ~x)) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ((~p) ^ 2, 1 / 4) && GtQ(~a, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * sin(~(e') + ~(f') * ~x) + ~(d') * cos(~(e') + ~(f') * ~x)) ^ 2) ^ ~p, ~x) => ((~a + ~b * (~c * sin(~e + ~f * ~x) + ~d * cos(~e + ~f * ~x)) ^ 2) ^ ~p / (1 + (~b * (~c * sin(~e + ~f * ~x) + ~d * cos(~e + ~f * ~x)) ^ 2) / ~a) ^ ~p) * Antiderivative((1 + (~b * (~c * sin(~e + ~f * ~x) + ~d * cos(~e + ~f * ~x)) ^ 2) / ~a) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ((~p) ^ 2, 1 / 4) && Not(GtQ(~a, 0)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.7 (d trig)^m (a+b (c sin)^n)^p.jl


text:
# This file was generated, do not modify it. # hide
describe(smarket, :mean, :std, :eltype)


text:
export solve

function possible(p::Puzzle, x::Int, y::Int, n::Int)::Bool
    for j in 1:9
        if p[x,j] == n
            return false;
        end
    end
    for i in 1:9
        if p[i,y] == n
            return false;
        end
    end
    # get sub-grid coordinates
    xs = Int.(3 * ceil(x / 3 - 1) .+ collect(1:3))
    ys = Int.(3 * ceil(y / 3 - 1) .+ collect(1:3))
    for i in xs
        for j in ys
            if p[i,j] == n
                return false;
            end
        end
    end
    return true;
end

function solve(p::Puzzle)::Bool
    for i in 1:9
        for j in 1:9
            if p[i,j] == 0
                for n in 1:9
                    if possible(p, i, j, n)
                        p[i,j] = n
                        # recursion
                        if solve(p)
                            return true;
                        else
                            # backtracking
                            p[i,j] = 0
                        end
                    end
                end
                return false;
            end
        end
    end
    return true;
end

max_stars_repo_name: Joshua-Robison/Sudoku

max_issues_repo_name: Joshua-Robison/Sudoku

max_forks_repo_name: Joshua-Robison/Sudoku


text:
"""Calculate the Fermi Function given some energy and some temperature.

Parameters
----------
en : float
    Energy relative to the fermi energy (E-Ef).

kbt : Temperature in same units as Energy (using kB).

Returns
-------
ffun : float
    The Fermi Function at en and temp."""
function fermi_fun(en, kbt)
    if en == 0
        ffun = 0.5
    elseif kbt == 0
        if en < 0
            ffun = 1.0
        elseif en > 0
            ffun = 0.0
        else
            throw(DomainError("Energy must be a real number."))
        end
    elseif kbt > 0
        #Using a tanh here instead of the more traditional formulation because
        #the function is better behaved for very large or small values of en/kbt
        ffun = 0.5*(1-tanh(0.5*en/kbt))

    elseif kbt < 0
        throw(DomainError("Temperature must be >= 0."))
    end

    return ffun
end


text:
module triangle_area

# package code goes here
export heron
@doc "
   헤론의 공식을 이용하여 세 변의 길이를 입력하면 삼각형의 넓이를 구해 주는 함수

   사용법 : n=heron(a,b,c)

   a,b,c : 삼각형 각 변의 길이

   n : 삼각형의 넓이"->
   function heron(a,b,c)
   s=(a+b+c)/2
   return sqrt(s*(s-a)*(s-b)*(s-c))
   end

end # module


text:
# The zoo society
abstract type Animal end
abstract type AbstractBird <: Animal end
abstract type AbstractCat <: Animal end

# We will define some required methods for Animal
# Some of them have a default returns

# All subtype of Animal should have a name
name(animal::Animal) = animal.name
gender(animal::Animal) = "unknown"

# We can know about each animal if it can fly
# By requiring this function (or property)
# Most animals won't fly
flyable(animal::Animal) = false
# Most birds can
flyable(bird::AbstractBird) = true

# play is a kind of relationship inside the society of animals
# Most animal won't play with others
play(a::Animal, b::Animal) = "$(name(a)) won't play with $(name(b))"
play(a::Animal, b::Animal, c::Animal) = "$(name(a)), $(name(b)), $(name(c)) won't play together"

# normal cats
struct Cat <: AbstractCat
    name::String
end

# normal birds
struct Bird <: AbstractBird
    name::String
end

# Kitty is a special Cat
struct Kitty <: AbstractCat
    weight::Float64
end

# Kitty is 0.3kg when she was born
Kitty() = Kitty(0.3)

# comment this method to see what will happen
# Kitty's name
name(me::Kitty) = "Kitty"
# Kitty's gender
gender(me::Kitty) = "female"

# Kohen is a special bird
struct Kohen <: AbstractBird
    ishome::Bool
end

# Kohen is usually not at home
Kohen() = Kohen(false)

# Kohen's name
name(me::Kohen) = "Kohen"
# Kohen's sex
gender(me::Kohen) = "male"

# Kohen cannot fly
flyable(me::Kohen) = false

# Kohen will play with Kitty
# comment this to see what will happen
play(bird::Kohen, cat::Kitty) = "$(name(bird)) will play with $(name(cat))"
play(cat::Kitty, bird::Kohen) = play(bird, cat)
# play(cat::Kitty, bird::Kohen) = "$(name(cat)) won't play with $(name(bird))"

# We will edit the print content of Kitty and Kohen
import Base: show

show(io::IO, kitty::Kitty) = print(io, "Kitty (weight $(kitty.weight))")
show(io::IO, kohen::Kohen) = print(io, "Kohen (is home: $(kohen.ishome))")

###################################
# Bear can have different colors
abstract type AbstractBearColor end
abstract type Green <: AbstractBearColor end
abstract type White <: AbstractBearColor end
abstract type Brown <: AbstractBearColor end
abstract type AbstractBear{T <: AbstractBearColor} <: Animal end

# we could define some general interface to All KINDS OF bears
# and use multiple dispatch to specify methods to certain
# type with certain type parameters (signature)
# the word where declares type parameters
color(::Type{AbstractBear{T}}) where T = "We don't have $T Bear"
# you can specify type parameters' parent type if you want
# Union will create a type union, behaves like an abstract type
# in type inference, but it cannot be subtyped
color(::Type{AbstractBear{T}}) where {T <: Union{Brown, White}} = T
color(bear::AbstractBear) where T = color(typeof(bear))

# Abram is a Brown bear
struct Abram <: AbstractBear{Brown}
end

name(bear::Abram) = "Abram"

# Dima is a White bear
struct Dima <: AbstractBear{White}
end

name(bear::Dima) = "Dima"

play(a::Abram, b::Dima, c::Kohen) = "Abram, Dima and Kohen often play together"
# * uncomment this to see what happen
# play(a::Abram, b::Dima, c::Kitty) = "Abram and Dima play with Kitty too"

# Playground

bird = Bird("bird")
cat = Cat("cat")
kohen = Kohen()
kitty = Kitty()
abram = Abram()
dima = Dima()

println("There is a cat and a bird in the Zoo")
println("The cat is ", kitty)
println("The bird is ", kohen)
println("Normally, ", play(bird, cat))
println("And, ", play(cat, bird))
println("However, ", play(kitty, kohen))
println("And, ", play(kohen, kitty))
print("Can $(name(kohen)) fly? ")
println(if flyable(kohen) "Yes." else "Nope." end)
println("That's why...")

println("The family of bear also lives in the Zoo")
println("Do we have Green Bear? ", color(AbstractBear{Green}))
println("What do we have? ", color(AbstractBear{Brown}), " Bear and ", color(AbstractBear{White}), " Bear.")
println(play(abram, dima, kohen))
println(play(abram, dima, kitty))


text:
# follow the steepest path to minima

struct EnergyDescent{T3} <: SpinSampler{T3}
    A::T3           # q × N × M
    A₁::T3          # 1 × N × M
    α::T3           # 1 × 1 × M
    q::Int          # number of spin states
    N::Int          # number of spins
    M::Int          # sample size
    function EnergyDescent(X::AbstractArray{<:AbstractFloat,3})
        q, N, M = size(X)
        _eltype = promote_type(eltype(X), Float32)
        A = similar(X, _eltype, q, N, M)
        A₁= similar(X, _eltype, 1, N, M)
        α = similar(X, _eltype, 1, 1, M)
        return new{typeof(A)}(A, A₁, α, q, N, M)
    end
end

(f::EnergyDescent)(X, Δenergy) = energydescent!(X, Δenergy, f)

function energydescent!(X::AbstractArray{T,3}, Δenergy, buffer = EnergyDescent(X)) where T
    (; A, A₁, α) = buffer
    _onehalf = convert(T, 1//2)

    ΔE = Δenergy(X)                 # E(xᵠᵏ)-E(x)
    δE, I = findmin(ΔE, dims = 1:2)

    fill!(A, zero(T))               # initialize
    fill!(@view(A[I]), one(T))      # onehot encode flipped spin and proposed spin state
    sum!(A₁, A)                     # onehot encode site of flipped spin
    @. A = X + (A - A₁ * X)         # proposed state
    @. A = A > _onehalf             # ensure the entries are in (0,1)

    @. α = (δE < 0) & (δE ≉ 0)      # flip to new state only if energy decreases

    @. X = (1-α) * X + α * A        # update X
    @. X = X > _onehalf             # ensure the entries are in (0,1)
    return X
end


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

"""
    PointSetData(data, coords)

Spatial `data` georeferenced with `coords`, which can be a matrix or
a vector of tuples. The `data` argument is a dictionary mapping variable
names to Julia arrays with the actual data.

See also: [`PointSet`](@ref)
"""
struct PointSetData{T,N} <: AbstractData{T,N}
  data::OrderedDict{Symbol,<:AbstractArray}
  domain::PointSet{T,N}

  function PointSetData{T,N}(data, domain) where {N,T}
    nvals = [length(array) for array in values(data)]
    @assert all(nvals .== npoints(domain)) "data and coords must have the same number of points"
    new(data, domain)
  end
end

PointSetData(data::OrderedDict{Symbol,<:AbstractArray},
             coords::AbstractMatrix{T}) where {T} =
  PointSetData{T,size(coords,1)}(data, PointSet(coords))

PointSetData(data::OrderedDict{Symbol,<:AbstractArray},
             coords::AbstractVector{NTuple{N,T}}) where {N,T} =
  PointSetData(data, [c[i] for i in 1:N, c in coords])

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::PointSetData{T,N}) where {N,T}
  npts = npoints(geodata)
  print(io, "$npts PointSetData{$T,$N}")
end


text:
doc"""
Spherically symmetric function

$f(r,\theta,\phi) = u(r)$

where $u(r)$ is a user-defined radial function.

The function vanishes for $r>r_{max}$.

f = RadialFunc(u; R=(0.,0.,0.), rmax=Inf)

`u` - a user-defined radial function

`R` - origin of spherical coordinate system
"""
type RadialFunc{T} <: ScalarFunc3d{Float64}
    u::T
    R::Point3{Float64}
    rmax::Float64
end

RadialFunc(u; R=(0.,0.,0.), rmax=Inf) = RadialFunc(u,R,rmax)

RadialFunc(r::AbstractVector, u::AbstractVector; R=(0.,0.,0.),
        rmax=Inf, npts::Int=100) =
    RadialFunc(Spline1D(r,u,npts),R,rmax)

typealias RadialFuncSpl RadialFunc{Spline1D{Float64,Vector{Float64}}}

function Base.writemime{T}(io::IO, mime::MIME"text/plain", f::RadialFunc{T})
    write(io, summary(f), "\n")
    write(io, "u(r) = ")
    writemime(io, mime, f.u)
    write(io, "\nR = ")
    writemime(io, mime, f.R)
    write(io, "\nrmax = ")
    writemime(io, mime, f.rmax)
end

function Base.call(f::RadialFunc, x::Float64, y::Float64, z::Float64)
    x0,y0,z0 = f.R
    r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2)
    r > f.rmax ? 0. : f.u(r)
end

Base.call(f::RadialFunc, p::Point3{Float64}) = f(p[1],p[2],p[3])
@vectorize2 Base.call RadialFunc Point3{Float64} Float64

function grad(f::RadialFunc{Spline1D}, x::Real, y::Real, z::Real)
    x0,y0,z0 = f.R
    r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2)
    dfdr = r > f.rmax ? 0. : deriv1(f.u,r)
    (dfdr,0.,0.)
end

grad(f::RadialFunc, p::Point3{Float64}) = grad(f,p[1],p[2],p[3])
@vectorize2 grad RadialFunc Point3{Float64} Point3{Float64}

∇ = grad

function laplacian(f::RadialFunc{Spline1D}, x::Real, y::Real, z::Real)
    x0,y0,z0 = f.R
    r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2)
    _,df,df2 = eval_deriv12(f.u,r)
    r > f.rmax ? 0. : 2*df/r + df2
end

laplacian(f::RadialFunc, p::Point3{Float64}) = laplacian(f,p[1],p[2],p[3])
@vectorize2 laplacian RadialFunc Point3{Float64} Float64

function setposition(f::RadialFunc, R::Point3{Float64})
    f.R = R
end

max_stars_repo_name: johncsnyder/Element


text:
using FFTW
using DifferentialEquations
using Plots

nu = 0.001 # Diffusion constant

# Define spatial domain
L = 20 # Length of domain
N = 1000 # Number of discretization points
dx = L/N
x = -L/2:dx:L/2-dx # Define x domain

# Define discrete wavenumbers
kappa = (2pi/L)*(-N/2:N/2-1)
kappa = fftshift(kappa) # Re-order fft wavenumbers

# Initial condition
u0 = sech.(x)

# Simulate PDE in spatial domain
f(u, p, t) = real.(-u.*ifft(im*kappa.*fft(u)) + nu*ifft(-(kappa.^2).*fft(u)))
prob = ODEProblem(f, u0, (0.0,2.5))
u = solve(prob)

# Plot solution in time
p1 = surface(hcat(u.u...))


text:
struct TransientMultiFieldCellField{A} <: TransientCellField
  cellfield::A
  derivatives::Tuple
  transient_single_fields::Vector{<:TransientCellField} # used to iterate
end

MultiFieldTypes = Union{MultiFieldCellField,MultiFieldFEFunction}

function TransientCellField(multi_field::MultiFieldTypes,derivatives::Tuple)
  transient_single_fields = _to_transient_single_fields(multi_field,derivatives)
  TransientMultiFieldCellField(multi_field,derivatives,transient_single_fields)
end

function get_data(f::TransientMultiFieldCellField)
  s = """
  Function get_data is not implemented for TransientMultiFieldCellField at this moment.
  You need to extract the individual fields and then evaluate them separatelly.

  If ever implement this, evaluating a `MultiFieldCellField` directly would provide,
  at each evaluation point, a tuple with the value of the different fields.
  """
  @notimplemented s
end

get_triangulation(f::TransientMultiFieldCellField) = get_triangulation(f.cellfield)
DomainStyle(::Type{TransientMultiFieldCellField{A}}) where A = DomainStyle(A)
num_fields(f::TransientMultiFieldCellField) = length(f.cellfield)
gradient(f::TransientMultiFieldCellField) = gradient(f.cellfield)
∇∇(f::TransientMultiFieldCellField) = ∇∇(f.cellfield)
change_domain(f::TransientMultiFieldCellField,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.cellfield,trian,target_domain)

# Get single index
function Base.getindex(f::TransientMultiFieldCellField,ifield::Integer)
  single_field = f.cellfield[ifield]
  single_derivatives = ()
  for ifield_derivatives in f.derivatives
    single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
  end
  TransientSingleFieldCellField(single_field,single_derivatives)
end

# Get multiple indices
function Base.getindex(f::TransientMultiFieldCellField,indices::Vector{<:Int})
  cellfield = MultiFieldCellField(f.cellfield[indices],DomainStyle(f.cellfield))
  derivatives = ()
  for derivative in f.derivatives
    derivatives = (derivatives...,MultiFieldCellField(derivative[indices],DomainStyle(derivative)))
  end
  transient_single_fields = _to_transient_single_fields(cellfield,derivatives)
  TransientMultiFieldCellField(cellfield,derivatives,transient_single_fields)
end

function _to_transient_single_fields(multi_field,derivatives)
  transient_single_fields = TransientCellField[]
  for ifield in 1:num_fields(multi_field)
    single_field = multi_field[ifield]
    single_derivatives = ()
    for ifield_derivatives in derivatives
      single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
    end
    transient_single_field = TransientSingleFieldCellField(single_field,single_derivatives)
    push!(transient_single_fields,transient_single_field)
  end
  transient_single_fields
end

# Iterate functions
Base.iterate(f::TransientMultiFieldCellField)  = iterate(f.transient_single_fields)
Base.iterate(f::TransientMultiFieldCellField,state)  = iterate(f.transient_single_fields,state)

# Time derivative
function ∂t(f::TransientMultiFieldCellField)
  cellfield, derivatives = first_and_tail(f.derivatives)
  transient_single_field_derivatives = TransientCellField[]
  for transient_single_field in f.transient_single_fields
    push!(transient_single_field_derivatives,∂t(transient_single_field))
  end
  TransientMultiFieldCellField(cellfield,derivatives,transient_single_field_derivatives)
end

∂tt(f::TransientMultiFieldCellField) = ∂t(∂t(f))


text:
using ReduceAlgebra
using Test, Pkg

# write your own tests here
@test (Pkg.test("Reduce"); true)
@test (Pkg.test("ReduceLinAlg"); true)


text:
module Wavefunctions

using LinearMaps
using RecipesBase

import Base: norm, normalize!, *, A_mul_B!

abstract type AbstractWfn{T} end

# Single-active electron wavefunctions

# An AbstractSAEWfn is assumed to have a data matrix, containing the
# wavefunction with increasing spatial dimension along the rows, and
# increasing partial wave index along the columns.
abstract type AbstractSAEWfn{T} <: AbstractWfn{T} end

function normalize!(wfn::AbstractSAEWfn)
    scale!(wfn.data, 1/norm(wfn))
    wfn
end

function A_mul_B!(Y::W, A::Union{AbstractMatrix,AbstractLinearMap}, B::W) where W <: AbstractSAEWfn
    A_mul_B!(view(Y.data, :), A, view(B.data, :))
    Y
end

include("wfn_utils.jl")
include("eigenstates.jl")
include("fd_sae_wfn.jl")

export AbstractWfn

end # module


text:
function crossprocess_wind(coil::CoilGeometry, Ω::Unitful.ElectricalResistance,
    processfill::Number,targetfill::Number)

    targetwind=optimalcoil(coil,Ω,targetfill)
    turns  =targetwind.turns
    gauge  = AWG_Chart[targetwind.gauge]
    targetΩ=targetwind.Ω

    f=id -> begin
        newcoil = CoilGeometry(id, 2coil.od,coil.len)
        estimatetruefill(newcoil,gauge;Ω=targetΩ,maxturns=turns,fillfactor=processfill).Ω-targetΩ
    end
    id=find_zero(f,coil.id)
    y=estimatetruefill(CoilGeometry(coil;id=id),gauge;Ω=targetΩ,maxturns=turns,fillfactor=processfill)
    merge((gauge=targetwind.gauge,),y)
end

"""
Calculate the minimum length coil necessary to achieve the minimum turns of wire

    minimumcoillength(coil, Ω, minturns, fillfactor=1.0)
    minimumcoillength(coil, v, p, minturns, fillfactor=1.0)

"""
function minimumcoillength(coil::CoilGeometry, Ω::Unitful.ElectricalResistance,
    minturns::Integer, fillfactor::Real=1.0)

    f=let Ω=u"Ω"(Ω)  #Simplify for type inference, necessary on <=1.3?
        x -> begin
            _coil = CoilGeometry(coil.id,coil.od,x)
            # _Ω = upreferred(Ω)
            optimalcoil(_coil, Ω, fillfactor)
        end
    end

    rng=5u"mm":0.01u"mm":300u"mm"
    m=mappedarray(x->f(x).turns,rng)
    i=searchsortedfirst(m,minturns)
    f(rng[i])
end

minimumcoillength(coil::CoilGeometry, v::Unitful.Voltage, p::Unitful.Power,
    minturns::Integer, fillfactor::Real=1.0) = minimumcoillength(coil, u"Ω"(v^2/p), minturns, fillfactor)

"""

    gaugebreakpoint(coil::CoilGeometry, Ω::Unitful.ElectricalResistance, len::Unitful.Length, fillfactor::Real=1.0)

Calculate length of coil with the nearest wire gauge break point, such that the geometry has the highest fill factor.

"""
function gaugebreakpoint(coil::CoilGeometry, Ω::Unitful.ElectricalResistance,
    len::Unitful.Length, fillfactor::Real=1.0)

    f=let Ω=u"Ω"(Ω)  #Simplify for type inference, necessary on <=1.3?
        x -> begin
            _coil = CoilGeometry(coil.id,coil.od,x)
            # _Ω = upreferred(Ω)
            optimalcoil(_coil, Ω, fillfactor)
        end
    end
    gauge = f(len).gauge

    rng=5u"mm":0.01u"mm":300u"mm"
    m=mappedarray(x->f(x).gauge,rng)
    a=rng[searchsortedfirst(m,gauge;rev=true)]
    b=rng[searchsortedlast(m,gauge;rev=true)]
    abs(len-a) < abs(len-b) ? a : b
end

function minimumturns(coil::CoilGeometry, power::Unitful.Power, fillfactor::Real=1.0; op=minimum,
    vref::Unitful.Voltage=24u"V", rng::AbstractVector{<:Unitful.Voltage} = 12.0u"V":1u"V":208u"V")

    tmap(rng) do v
        Ω = v^2/power |> u"Ω"
        y = optimalcoil(coil, Ω, fillfactor)
        # @info (v=v,turns=y.turns,relative = y.turns * vref/v)
        y.turns * vref/v |> floor |> Int # Equivalent number of turns
    end |> op

    # g = v -> begin
    #     Ω = v^2/power |> u"Ω"
    #     y = optimalcoil(coil, Ω, fillfactor)
    #     # @info (v=v,turns=y.turns,relative = y.turns * vref/v, id=Threads.threadid())
    #     y.turns * vref/v |> floor |> Int
    # end
    # reduce(op, Map(g), rng)
end

"""

Calculate the coil wire length and weight from an `optimalcoil` calculation.

"""
function weight(a)
    coil=a.coil
    turns=a.turns
    Ω=a.Ω
    gauge=AWG_Chart[a.gauge]
    feet=Ω/gauge.Ω
    return (length=feet,weight=feet/gauge.ft_lb)
end


text:
module ASP

using QRupdate

export bpdual, as_bpdn, as_nnls, OptimResults

include("bpdual.jl")
include("solvers.jl")

end


text:
export time_fractional_op, time_fractional_t_op

@doc raw"""
    time_fractional_t_op(i::Union{Integer, PyObject}, ta::PyObject,
        α::Union{Float64, PyObject}, Δt::Union{Float64, PyObject})

Returns the coefficients for the time fractional derivative
```math
{}_0^CD_t^\alpha f(t) = \frac{1}{\Gamma(1-\alpha)}\int_0^t \frac{f'(\tau)d\tau}{(t-\tau)^\alpha}
```

The discretization scheme used here is

$\begin{aligned}
& {}_0^CD_\tau^\alpha u(\tau_n) \\
= &\frac{\Delta \tau^{-\alpha}}{\Gamma(2-\alpha)}\left[G_0 u_n - \sum_{k=1}^{n-1}(G_{n-k-1}-G_{n-k})u_k + G_n u_0 \right] + \mathcal{O}(\Delta \tau^{2-\alpha})
\end{aligned}$

Here
```math
G_m = (m+1)^{1-\alpha} - m^{1-\alpha}, \quad m\geq 0, \quad 0<\alpha<1
```

The function returns
- `c` : $\frac{\Delta \tau^{-\alpha}}{\Gamma(2-\alpha)}$
- `cum` : $- \sum_{k=1}^{n-1}(G_{n-k-1}-G_{n-k})u_k + G_n u_0$
"""
function time_fractional_t_op(i::Union{Integer, PyObject}, ta::PyObject,
         α::Union{Float64, PyObject}, Δt::Union{Float64, PyObject})
    α = convert_to_tensor(α, dtype=Float64)
    i = convert_to_tensor(i, dtype=Int32)

    function aki(k,i)
        a = cast(i-k, Float64)
        return (a+2.0)^(1-α)-2*(a+1.0)^(1-α)+a^(1-α)
    end

    function cond1(k, i, ta, cum)
        k<=i
    end

    function body1(k, i, ta, cum)
        u = read(ta, k)
        return k+1, i, ta, cum+aki(k,i)*u
    end

    cum = tf.zeros_like(read(ta, 1))
    k = constant(1, dtype=Int32)
    _, _, _, cum = while_loop(cond1, body1, [k, i, ta, cum])
    c = Δt^(-α)/exp(tf.math.lgamma(2-α))
    return c, cum
end

@doc raw"""
    time_fractional_op(α::Union{Float64, PyObject}, f_fun::Function,
        T::Float64, u0::Union{Float64, Array{Float64}, PyObject},
        NT::Int64, θ=missing)

Returns a $(NT+1)\times \texttt{size}(u_0)$ solution array. The function solves the following time fractional differential equation with explicit scheme

```math
{}_0^CD_t^\alpha u(t) = f(t, u, \theta)
```
"""
function time_fractional_op(α::Union{Float64, PyObject}, f_fun::Function,
        T::Float64, u0::Union{Float64, Array{Float64}, PyObject},
        NT::Int64, θ=missing)
        Δt = T/NT
        ta = TensorArray(NT+1)
        ta = write(ta, 1, convert_to_tensor(u0))
        function condition(i, ta)
            i<=NT+1
        end
        function body(i, ta)
            # time = (i-1)*Δt,
            u = read(ta, i-1)
            k = cast(i,Float64)-1
            c, cum = time_fractional_t_op(k, ta, α, Δt)
            F = f_fun(k*Δt, u, θ)
            ta = write(ta, i, F/c - cum)
            i+1, ta
        end
        i = constant(2, dtype=Int32)
        _, out = while_loop(condition, body, [i, ta])
        return stack(out)
end