text:
struct QuadraticMeltRate{T <: Real} <: AbstractMeltRate
    γT :: T                       #(calibrated) heat exchange velocity
    λ1 :: T                       #liquidus slope
    λ2 :: T                       #liquidus intercept
    λ3 :: T                       #liquidus pressure coefficient
    ρi :: T                       #ice density
    ρw :: T                       #water density
    L :: T                        #latent heat of fusion for ice
    c :: T                        #specific heat capacity of ocean water
    Ta :: Function                #Function specifying ambient temperature
    Sa :: Function                #Function specifying ambient salinity
    flocal :: Bool                #Flag for local or non-local temperature dependence
    melt_partial_cell :: Bool       #Flag for melting applied to partial cells or not
end

"""
    QuadraticMeltRate(; ,kwargs)

Construct a QuadraticMeltRate object to prescribe the melt rate in WAVI

Keyword arguments
=================
- γT: calibrated heat exchange velocity (units m/s).
- λ1: liquidus slope (units ∘C)
- λ2: liquidus intercept (units ∘C)
- λ3: liquidus pressure coefficient (units K/m)
- ρi: ice density (units kg/m^3)
- ρw: water density (units kg/m^3)
- L: latent heat of fusion for ice (units: J/kg)
- c: specific heat capacity of water (units J/kg/K)
- Ta: ambient temperature function, defaults to the ISOMIP warm0 scenario. Ta must be a function of a single variable (depth) [time dependence not yet implemented in WAVI.jl]
- Sa: ambient salnity function, defaults to the ISOMIP warm0 scenario.  Sa must be a function of a single variable (depth) [time dependence not yet implemented in WAVI.jl]
- flocal: flag for local or non-local temperature dependence. Set to true for local dependence (eqn 4 in Favier 2019 GMD) or false for non-local dependence (eqn 5)
- melt_partial_cell: flag for specify whether to apply melt on partially floating cells (true) or not (false)
"""

function QuadraticMeltRate(;
                            γT = 1.e-3,
                            λ1 = -5.73e-2,
                            λ2 = 8.32e-4,
                            λ3 = 7.61e-4,
                            L = 3.35e5,
                            c = 3.974e3,
                            ρi = 918.8,
                            ρw = 1028.0,
                            Ta = isomip_warm0_temp,
                            Sa = isomip_warm0_salinity,
                            flocal = true,
                            melt_partial_cell = false)

    #check that Ta and Sa accept a single argument
    try Ta(0.0)
    catch
        ArgumentError("Ambient temperature function Ta must accept a single argument")
    end
    try Sa(0.0)
    catch
        ArgumentError("Ambient salinity function Sa must accept a single argument")
    end

    return QuadraticMeltRate(γT, λ1, λ2, λ3, ρi, ρw, L, c, Ta, Sa, flocal, melt_partial_cell)
end

"""
    update_melt_rate!(quad_melt_rate::QuadraticMeltRate, fields, grid)

Wrapper script to update the melt rate for a QuadraticMeltRate.
"""
function update_melt_rate!(quad_melt_rate::QuadraticMeltRate, fields, grid)
    @unpack basal_melt, h, b, grounded_fraction = fields.gh #get the ice thickness and grounded fraction

    #compute the ice draft
    zb = fields.gh.b .* (grounded_fraction .== 1) + - quad_melt_rate.ρi / quad_melt_rate.ρw .* fields.gh.h .* (grounded_fraction .< 1)

    set_quadratic_melt_rate!(basal_melt,
                            quad_melt_rate,
                            zb,
                            grounded_fraction)
    return nothing
end

function set_quadratic_melt_rate!(basal_melt,
                                qmr,
                                zb,
                                grounded_fraction)

    #compute freezing temperature and Ta - Tf everywhere
    Sa_shelf = qmr.Sa.(zb)
    Ta_shelf = qmr.Ta.(zb)
    Tstar = qmr.λ1 .* Sa_shelf .+ qmr.λ2 .+ qmr.λ3 .* zb .- Ta_shelf
    Tstar_shelf_mean = sum(Tstar[grounded_fraction .== 0])/length(Tstar[grounded_fraction .== 0])

    #set melt rate
    if (qmr.melt_partial_cell) && (qmr.flocal) #partial cell melting and local
        basal_melt[:] .=  qmr.γT .* (qmr.ρw * qmr.c / qmr.ρi /qmr.L)^2 .* Tstar[:].^2 .* (1 .- grounded_fraction[:])

    elseif ~(qmr.melt_partial_cell) && (qmr.flocal) #no partial cell melting and local
        basal_melt[grounded_fraction .== 0] .=  qmr.γT .* (qmr.ρw * qmr.c / qmr.ρi /qmr.L)^2 .* Tstar[grounded_fraction .== 0].^2
        basal_melt[.~(grounded_fraction .== 0)] .= 0

    elseif (qmr.melt_partial_cell) && ~(qmr.flocal) #partial cell melting and non-local
        basal_melt[:] .=  qmr.γT .* (qmr.ρw * qmr.c / qmr.ρi /qmr.L)^2 .* Tstar[:].^2 .* (1 .- grounded_fraction[:])

    elseif ~(qmr.melt_partial_cell) && ~(qmr.flocal) #no partial cell and non-local
        basal_melt[grounded_fraction .== 0] .=  qmr.γT .* (qmr.ρw * qmr.c / qmr.ρi /qmr.L)^2 .* Tstar[grounded_fraction .== 0] .* Tstar_shelf_mean
        basal_melt[.~(grounded_fraction .== 0)] .= 0
    end
    basal_melt[:] .= basal_melt[:].* 365.25*24*60*60
    return nothing
end


text:
function isemsa(io::IO)
	for (lx, line) in enumerate(eachline(io))
		tmp = split_emsa_header_item(line)
		if isnothing(tmp) || #
		   (lx==1 && !(tmp[1]=="FORMAT" && uppercase(tmp[2])=="EMSA/MAS SPECTRAL DATA FILE")) || #
		   (lx==2 && !(tmp[1]=="VERSION" && tmp[2]=="1.0"))
			return false
		end
		if lx>=2
			return true
		end
	end
end

isemsa(filename::AbstractString) =
	open(filename,"r") do f
		isemsa(f)
	end

function split_emsa_header_item(line::AbstractString)
    p = findfirst(":",line)
    if p ≠ nothing
        tmp = uppercase(strip(SubString(line,2:p.start-1)))
        pp = findfirst("-",tmp)
        if pp ≠ nothing
            key, mod = strip(SubString(tmp,1:pp.start-1)), strip(SubString(tmp,pp.stop+1))
        else
            key, mod = tmp, nothing
        end
        value = strip(SubString(line,p.stop+1))
        # println(key," = ",value)
        return (key, value, mod)
    else
        return nothing
    end
end

function parsecoating(value::AbstractString)::Film
    i=findfirst(" nm of ",value)
	if !isnothing(i)
		try
			thk = parse(Float64, value[1:i.start-1])*1.0e-7 # cm
			mat = parsedtsa2comp(value[i.stop+1:end])
			return Film(mat, thk)
		catch err
			@warn "Error parsing $(value) as a coating."
		end
	end
    return Film()
end

"""
    readEMSA(filename::AbstractString, T::Type{<:Real}=Float64)::Spectrum

Read an ISO/EMSA format spectrum from a disk file at the specified path.
T is the type of the channel data elements.
"""
function readEMSA(filename::AbstractString, T::Type{<:Real}=Float64)::Spectrum
	open(filename,"r") do f
		res = readEMSA(f, T)
		_setfilename(res, filename)
		return res
	end
end

function _setfilename(spec::Spectrum, filename::String)
	spec[:Filename]=filename
	if startswith( spec[:Name], "Bruker")
		if haskey(spec, :Composition)
			spec[:Name]=name(spec[:Composition])
		else
			path = splitpath(filename)
			spec[:Name]=path[end]
		end
	end
	if endswith(spec[:Name],".txt") || endswith(spec[:Name],".msa")
		spec[:Name] = spec[:Name][1:end-4]
	end
end

"""
    readEMSA(ios::IO, T::Type{<:Real}=Float64)::Spectrum

Read an ISO/EMSA format spectrum from a disk file at the specified path.
T is the type of the channel data elements.
"""
function readEMSA(f::IO, T::Type{<:Real}=Float64)::Spectrum
    energy, counts, eaxis = LinearEnergyScale(0.0,10.0), T[], Float64[]
    props = Dict{Symbol,Any}()
    props[:Filename]=filename
    inData, xpcscale = 0, 1.0
    stgpos = Dict{Symbol,Float64}()
	date, time, datatype = missing, missing, :Y
    for (lx, line) in enumerate(eachline(f))
        if (lx ≤ 2)
            res = split_emsa_header_item(line)
            if isnothing(res) || ((lx==1) && (res[1]!="FORMAT" || uppercase(res[2])!="EMSA/MAS SPECTRAL DATA FILE"))
                error("This file does not have the correct header to be an EMSA/MAS spectrum file.")
            end
            if isnothing(res) || ((lx==2) && (res[1]!="VERSION" || res[2]!="1.0"))
                error("This file is not a VERSION=1.0 EMSA/MAS spectrum file.")
            end
        elseif inData>0
            if startswith(line, "#ENDOFDATA")
                break
            else
				if datatype == :Y
					for cc in filter(s->length(s)>0, strip.(split(line, ",")))
						try
							append!(counts, parse(T, cc))
						catch
							@warn "Unparsable counts while reading in \"$cc\" from \"$line\"."
						end
						inData+=1
					end
				else
					@assert datatype==:XY
					try
						items = filter(s->length(s)>0, strip.(split(line, ",")))
						for i in 1:2:length(items)
							append!(eaxis, parse(Float64, items[i]))
							append!(counts, parse(T, items[i+1]))
							inData+=1
						end
					catch
						@warn "Unparsable counts while reading \"$line\"."
					end

				end
			end
        elseif startswith(line,"#")
			try
	            res = split_emsa_header_item(line)
	            if res ≠ nothing
	                key, value, mod = res
	                if key == "SPECTRUM"
	                    inData = 1
	                elseif key == "BEAMKV"
	                    props[:BeamEnergy]=1000.0*parse(Float64,value) # in eV
	                elseif key == "XPERCHAN"
						xperch=parse(Float64,value)
						xpcscale = isnothing(mod) ? (xperch < 0.1 ? 1000.0 : 1.0) :  # Guess likely eV or keV
										(isequal(mod,"KEV") ? 1000.0 : 1.0)
	                    energy = LinearEnergyScale(xpcscale*energy.offset, xpcscale*xperch)
	                elseif key == "LIVETIME"
	                    props[:LiveTime]=parse(Float64,value)
	                elseif key == "REALTIME"
	                    props[:RealTime]=parse(Float64,value)
	                elseif key == "PROBECUR"
	                    props[:ProbeCurrent]=parse(Float64,value)
	                elseif key == "OFFSET"
	                    energy = LinearEnergyScale(xpcscale*parse(Float64,value), energy.width)
					elseif key == "DATE"
						try
							date = Date(parse(DateTime, value, DateFormat("d-u-y","english")))
						catch
							date = Date(parse(DateTime, value, DateFormat("d-m-y","english")))
						end
					elseif key == "TIME"
						time = Time(DateTime(value, count(c->c==':',value) == 2 ? "H:M:S" : "H:M"))
	                elseif key == "TITLE"
						props[:Name]=value
						props[:Title]=value
	                elseif key == "OWNER"
	                    props[:Owner]=value
	                elseif key == "ELEVANGLE"
	                    props[:TakeOffAngle] = (props[:Elevation] = deg2rad(parse(Float64, value)))
	                elseif key == "XPOSITION"
	                    stgpos[:X] = parse(Float64, value)
	                elseif key == "YPOSITION"
	                    stgpos[:Y] = parse(Float64, value)
	                elseif key == "ZPOSITION"
	                    stgpos[:Z] = parse(Float64, value)
	                elseif key == "COMMENT"
	                    prev = getindex(props,:Comment,missing)
	                    props[:Comment] = prev ≠ missing ? prev*"\n"*value : value
					elseif key == "DATATYPE"
						datatype = uppercase(value)=="XY" ? :XY : :Y
	                elseif key == "#D2STDCMP"
	                    props[:Composition] = parsedtsa2comp(value)
	                elseif key == "#CONDCOATING"
	                    props[:Coating] = parsecoating(value)
	                elseif key =="#RTDET" # Bruker
						props[:DetectorModel] = value
					elseif key == "#TDEADLYR" # Bruker
						@assert isnothing(mod) || (!isequal(mod,"CM")) "Unexpected scale -$mod"
						props[:DeadLayer] = parse(Float64, value)
					elseif key == "#FANO" # Bruker
						props[:Fano] = parse(Float64, value)
					elseif key == "#WORKING"
						props[:WorkingDistance] = 0.1*parse(Float64, split(value)[1])
					elseif key == "#MNFWHM" # Bruker
						sc = (isnothing(mod) ? 1000.0 : (isequal(mod,"KEV") ? 1000.0 : 1.0))
						props[:FWHMMnKa] = sc*parse(Float64, value)
					elseif key == "#IDENT" # Bruker
						elm = nothing
						try
							props[:XRFAnode] = parse(Element, value)
						catch
							# Just ignore it
						end
					elseif key== "#OXINSTELEMS"
						try
							props[:Elements] = map(v->parse(Element, v), strip.(split(value,",")))
						catch
							# Just ignore it
						end
					end
	            end
			catch
				println(stderr,"Error parsing: $line")
			end
        end
    end
	if !(ismissing(date) || ismissing(time))
		props[:AcquisitionTime] = DateTime(date, time)
	end
    if startswith( props[:Name], "Bruker") && haskey(props, :Composition)
        props[:Name]=name(props[:Composition])
    end
    if !isempty(stgpos)
		props[:StagePosition] = stgpos
	end
    return Spectrum(energy, counts, props)
end

"""
    readEMSA(filename::AbstractString, det::Detector, force::Bool=false)::Spectrum

Read an EMSA file and apply the specified detector.  If force is false and the detector and
read calibration don't match then the function errors.
"""
function readEMSA(filename::AbstractString, det::Detector, force::Bool=false, T::Type{<:Real}=Float64)::Spectrum
	spec = readEMSA(filename, T)
	if !force
		if abs(channel(energy(1, det),spec)-1)>1
			error("Spectrum and detector calibrations don't match - Offset.")
		end
		test = det.channelcount÷2
		if abs(channel(energy(test,det),spec)-test)>1
			error("Spectrum and detector calibrations don't match - Gain.")
		end
	end
	return apply(spec, det)
end

function writeEMSA(filename::AbstractString, spec::Spectrum)
	open(filename,"w") do ios
		writeEMSA(ios, spec)
	end
end

function writeEMSA(io::IOStream, spec::Spectrum)
	forceascii(ss) = map(c-> Int(c)&0x80==0 ? c : '?', ss)
	# "#FORMAT      : EMSA/MAS Spectral Data File"
	writeline(io, tag, data, extra="") =
		println(io,"#$(tag)$(repeat(' ',12-(length(tag)+length(extra))))$(extra): $(forceascii(data))")
	writeline(io, "FORMAT","EMSA/MAS Spectral Data File")
	writeline(io, "VERSION","1.0")
	writeline(io, "TITLE",spec[:Name])
	if haskey(spec,:AcquisitionTime)
		dt = spec[:AcquisitionTime]
		writeline(io, "DATE",uppercase(Dates.format(dt,"d-u-yyyy",locale="english"))) # 25-FEB-2022
		writeline(io, "TIME",Dates.format(dt,"HH:MM:SS")) # 15:32:23
	end
	if haskey(spec,:Owner)
		writeline(io, "OWNER",spec[:Owner])
	end
	writeline(io, "NPOINTS","$(length(spec))")
	writeline(io, "NCOLUMNS","1")
	writeline(io, "XUNITS","eV")
	writeline(io, "YUNITS","counts")
	if spec.energy isa LinearEnergyScale
		writeline(io, "DATATYPE","Y")
		writeline(io, "XPERCHAN","$(spec.energy.width)")
		writeline(io, "OFFSET", "$(spec.energy.offset)")
	else
		writeline(io, "DATATYPE","XY")
		dx = (energy(length(spec),spec)-energy(1,spec))/(length(spec)-1.0) # Mean width
		writeline(io, "XPERCHAN","$(dx)")
		writeline(io, "OFFSET", "$(energy(spec,1))") # energy of first channel
	end
	writeline(io, "SIGNALTYPE","EDS")
	if haskey(spec,:LiveTime)
		writeline(io, "LIVETIME","$(spec[:LiveTime])")
	end
	if haskey(spec,:RealTime)
		writeline(io, "REALTIME","$(spec[:RealTime])")
	end
	if haskey(spec,:BeamEnergy)
		writeline(io, "BEAMKV","$(0.001*spec[:BeamEnergy])")
	end
	if haskey(spec,:ProbeCurrent)
		writeline(io,"PROBECUR","$(spec[:ProbeCurrent])")
	end
	if haskey(spec,:TakeOffAngle)
		writeline(io, "ELEVANGLE","$(rad2deg(spec[:TakeOffAngle]))")
	end
	if haskey(spec,:Azimuthal)
		writeline(io, "AZIMANGLE","$(rad2deg(spec[:Azimuthal]))")
	end
	if haskey(spec,:StagePosition)
		sp=spec[:StagePosition]
		for (tag, sym) in ( ("XPOSITION", :X ), ("YPOSITION", :Y ), ("ZPOSITION", :Z ) )
			if haskey(sp, sym)
				writeline(io, tag, "$(10.0*sp[sym])","mm")
			end
		end
	end
	if haskey(spec,:Composition)
		writeline(io, "#D2STDCMP",todtsa2comp(spec[:Composition]))
	end
	if haskey(spec,:Specimen)
		writeline(io,"#SPECIMEN","$(spec[:Specimen])")
	end
	if haskey(spec,:Coating)
		cc = props[:Coating]
		writeline(io, "#CONDCOATING", "$(1.0e7*cc.thickness) nm of $(todtsa2comp(cc.material))")
	end
	if haskey(spec,:DetectorModel)
		writeline(io, "#RTDET",spec[:DetectorModel]) # Bruker
	end
	if haskey(spec,:DeadLayer)
		writeline(io, "#TDEADLYR", spec[:DeadLayer])  # Bruker
	end
	if haskey(spec,:Fano)
		writeline(io, "#FANO",spec[:Fano]) # Bruker
	end
	if haskey(spec, :FWHMMnKa)
		writeline(io, "#MNFWHM", spec[:FWHMMnKa], "EV") # Bruker
	end
	if haskey(spec, :XRFAnode)
		writeline(io, "#IDENT", spec[:XRFAnode]) # Bruker
	end
	if spec.energy isa LinearEnergyScale
		writeline(io, "SPECTRUM","Spectrum data follows in counts")
		for i in eachindex(spec)
			println(io, spec.counts[i], ",")
		end
	else
		writeline(io, "SPECTRUM","Spectrum data follows in energy, counts")
		for i in eachindex(spec)
			println(io, energy(i,spec), ", ", spec.counts[i], ",")
		end
	end
	writeline(io, "ENDOFDATA","")
end

ext: jl

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using PackageCompiler
using GaussianBeamLens
gui_file = joinpath(dirname(pathof(GaussianBeamLens)),"GtkGUI.jl")
#(@v1.6) pkg> activate .
#  Activating environment at `C:\Users\kittisopikulm\Documents\Julia\GaussianBeamLens\Project.toml`
create_sysimage([:MappedArrays, :IndirectArrays, :ColorSchemes, :Images, :GaussianBeamLens, :GtkReactive, :Gtk, :ImageView, :FFTW], precompile_execution_file=gui_file, sysimage_path="gaussian_beam_lens_gui.dll")
#PS C:\Users\kittisopikulm\Documents\Julia\GaussianBeamLens> julia --sysimage .\gaussian_beam_lens_gui.dll --project=. .\src\GtkGUI.jl
#  1.318605 seconds (1.82 M allocations: 168.221 MiB, 13.45% gc time, 69.55% compilation time)

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
using Printf

include("alias.jl")
include("utils.jl")
include("PGBF.jl")
include("CGBF.jl")
include("overlap.jl")
include("kinetic.jl")

function test_kinetic()

  s = PGBF(1.0)
  c = init_CGBF( 0.0, 0.0, 0.0 )
  push!(c,1.0,1.0)

  @assert isapprox(kinetic(1.0, (0.0,0.0,0.0), (0,0,0),
                           1.0, (0.0,0.0,0.0), (0,0,0)), 2.9530518648229536)
  @assert isapprox( kinetic(s,s), 1.5 )
  #@assert isapprox(kinetic(c,c),1.5)
  @printf "test_kinetic is passed\n"
end

test_kinetic()


text:
############################################
# Julia 1.2 temporary fix - Julia PR 33303 #
############################################
if VERSION == v"1.2"
    @eval function namedtuple(::Type{NamedTuple{names, T}}, args::Tuple) where {names, T <: Tuple}
        if length(args) != length(names)
            throw(ArgumentError("Wrong number of arguments to named tuple constructor."))
        end
        # Note T(args) might not return something of type T; e.g.
        # Tuple{Type{Float64}}((Float64,)) returns a Tuple{DataType}
        $(Expr(:splatnew, :(NamedTuple{names,T}), :(T(args))))
    end
    @generated function ntmerge(nt1::NamedTuple{names1, T1}, nt2::NamedTuple{names2, T2}) where {names1, T1 <: Tuple, names2, T2 <: Tuple}
        names = (names1..., names2...)
        T = Tuple{T1.types..., T2.types...}
        f = :(NamedTuple{$names, $T})
        args = :((Tuple(nt1)..., Tuple(nt2)...))
        quote
            $(Expr(:splatnew, f, args))
        end
    end
else
    function namedtuple(::Type{NamedTuple{names, T}}, args::Tuple) where {names, T <: Tuple}
        return NamedTuple{names, T}(args)
    end
    ntmerge(nt1::NamedTuple, nt2::NamedTuple) = merge(nt1, nt2)
end

#####################################################
# Helper functions for vectorize/reconstruct values #
#####################################################

vectorize(d::UnivariateDistribution, r::Real) = [r]
vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = copy(r)
vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = copy(vec(r))

# NOTE:
# We cannot use reconstruct{T} because val is always Vector{Real} then T will be Real.
# However here we would like the result to be specifric type, e.g. Array{Dual{4,Float64}, 2},
# otherwise we will have error for MatrixDistribution.
# Note this is not the case for MultivariateDistribution so I guess this might be lack of
# support for some types related to matrices (like PDMat).
reconstruct(d::UnivariateDistribution, val::AbstractVector) = val[1]
reconstruct(d::MultivariateDistribution, val::AbstractVector) = copy(val)
function reconstruct(d::MatrixDistribution, val::AbstractVector)
    return reshape(copy(val), size(d))
end
function reconstruct!(r, d::Distribution, val::AbstractVector)
    return reconstruct!(r, d, val)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector)
    r .= val
    return r
end
function reconstruct(d::Distribution, val::AbstractVector, n::Int)
    return reconstruct(d, val, n)
end
function reconstruct(d::UnivariateDistribution, val::AbstractVector, n::Int)
    return copy(val)
end
function reconstruct(d::MultivariateDistribution, val::AbstractVector, n::Int)
    return copy(reshape(val, size(d)[1], n))
end
function reconstruct(d::MatrixDistribution, val::AbstractVector, n::Int)
    tmp = reshape(val, size(d)[1], size(d)[2], n)
    orig = [tmp[:, :, i] for i in 1:size(tmp, 3)]
    return orig
end
function reconstruct!(r, d::Distribution, val::AbstractVector, n::Int)
    return reconstruct!(r, d, val, n)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector, n::Int)
    r .= val
    return r
end


text:
"""
Primary driver for network pairwise.
"""
function network_pairwise(T, V, cfg)::Matrix{T}

    # Get input
    networkdata = get_network_data(T, V, cfg)

    # Get compute flags
    flags = get_network_flags(cfg)

    # Compute graph data based on compute flags
    graphdata = compute_graph_data(networkdata, cfg)

    # Send to main kernel
    single_ground_all_pairs(graphdata, flags, cfg)
end

function compute_graph_data(data::NetworkData{T,V}, cfg)::GraphData{T,V} where {T,V}

    i,j,v = data.coords

    idx = findfirst(x -> x < 1, i)
    idx != nothing && throw("Indices no good")
    idx = findfirst(x -> x < 1, j)
    idx != nothing && throw("Indices no good")

    m = max(i[end], j[end])
    A = sparse(i,j,v,m,m)
    A = A + A'

    cc = connected_components(SimpleWeightedGraph(A))

    t = @elapsed G = laplacian(A)
    csinfo("Time taken to construct graph laplacian = $t", cfg["suppress_messages"] in TRUELIST)

    # T = eltype(i)
    exclude_pairs = Tuple{V,V}[]

    nodemap = Matrix{V}(undef,0,0)
    polymap = Matrix{V}(undef,0,0)
    hbmeta = RasterMeta()
    cellmap = Matrix{T}(undef,0,0)

    cum = initialize_cum_vectors(v)

    GraphData(G, cc, data.fp, data.fp,
                exclude_pairs, nodemap, polymap, hbmeta, cellmap, cum)
end

function get_network_flags(cfg)

    # Computation flags
    is_raster = false
    is_advanced = cfg["scenario"] in ADVANCED
    is_alltoone = false
    is_onetoall = false
    grnd_file_is_res = cfg["ground_file_is_resistances"] in TRUELIST
    policy = Symbol(cfg["remove_src_or_gnd"])
    solver = cfg["solver"]

    # Output flags
    write_volt_maps = cfg["write_volt_maps"] in TRUELIST
    write_cur_maps = cfg["write_cur_maps"] in TRUELIST
    write_cum_cur_maps_only = cfg["write_cum_cur_map_only"] in TRUELIST
    write_max_cur_maps = cfg["write_max_cur_maps"] in TRUELIST
    set_null_currents_to_nodata = cfg["set_null_currents_to_nodata"] in TRUELIST
    set_null_voltages_to_nodata = cfg["set_null_voltages_to_nodata"] in TRUELIST
    compress_grids = cfg["compress_grids"] in TRUELIST
    log_transform_maps = cfg["log_transform_maps"] in TRUELIST

    o = OutputFlags(write_volt_maps, write_cur_maps,
                    write_cum_cur_maps_only, write_max_cur_maps,
                    set_null_currents_to_nodata, set_null_voltages_to_nodata,
                    compress_grids, log_transform_maps)

    NetworkFlags(is_raster, is_advanced, is_alltoone, is_onetoall,
                grnd_file_is_res, policy, solver, o)
end

struct NetworkFlags
    is_raster::Bool
    is_advanced::Bool
    is_alltoone::Bool
    is_onetoall::Bool
    grnd_file_is_res::Bool
    policy::Symbol
    solver::String
    outputflags::OutputFlags
end


text:
using KNITRO, JuMP, Base.Test

for m in [Model(solver=KnitroSolver()),
          Model(solver=KnitroSolver(KTR_PARAM_ALG=5)),
          Model(solver=KnitroSolver(hessopt=1)),
          Model(solver=KnitroSolver(options_file=joinpath(dirname(@__FILE__),
                                                          "tuner-fixed.opt"))),
          Model(solver=KnitroSolver(tuner_file=joinpath(dirname(@__FILE__),
                                                        "tuner-explore.opt")))]
    @variable(m, x[1:3]>=0)
    @NLobjective(m, Min, 9.0 - 8.0*x[1] - 6.0*x[2] - 4.0*x[3]
                            + 2.0*x[1]^2 + 2.0*x[2]^2 + x[3]^2
                            + 2.0*x[1]*x[2] + 2.0*x[1]*x[3])
    @constraint(m, x[1] + x[2] + 2.0*x[3] <= 3)
    solve(m)
    @test_approx_eq_eps getvalue(x[1]) 1.3333333 1e-5
    @test_approx_eq_eps getvalue(x[2]) 0.7777777 1e-5
    @test_approx_eq_eps getvalue(x[3]) 0.4444444 1e-5
    @test_approx_eq_eps getobjectivevalue(m) 0.1111111 1e-5
end


text:
using Base.Test
using ExclusionProcesses

@testset "sep" begin

# A SEP update either does nothing or changes two spins
spin_state = SingleSpinState(10,5)
init_spins = copy(spin_state.spins)
update!(SEP, OpenBC, spin_state)
diff = xor.(init_spins, spin_state.spins)
@test sum(diff) == 2 || sum(diff) == 0

# For open boundary conditions the changed spins must be adjacent
@test diff[indmax(diff) + 1] == true || sum(diff) == 0

# Update of multispin configuration must cause convergence.
# This means that each `UInt64` must change in a way that increases the absolute
# value of the deviation from 32 of the number of 1s in its binary representation.
score(ms::UInt64) = abs(count_ones(ms) - 32)
ms_state = MultiSpinState(10,5)
init_scores = [score(site) for site in ms_state.spins]
update!(SEP, OpenBC, ms_state)
final_scores = [score(site) for site in ms_state.spins]
@test all(final_scores .>= init_scores)

end # testset

max_stars_repo_name: AustenLamacraft/ExclusionProcesses

max_issues_repo_name: AustenLamacraft/ExclusionProcesses

max_forks_repo_name: AustenLamacraft/ExclusionProcesses


text:
module TokyoVariantAnalysis4B2
using Logging
using DataFrames
using Dates

using ..TokyoVariantAnalysis4: date_to_value, value_to_date
using ..TokyoVariantAnalysis4A

export heatmap_array
export hm_x_min, hm_x_max, hm_y_min, hm_y_max, hm_x_size, hm_y_size, hm_f

@info "module", @__MODULE__

@info "--------"
@info "ヒートマップ形式尤度設定"
const hm_x_size = 500
const hm_y_size = 500
const hm_x_min = date_to_value(Date("2021-05-25"))
const hm_x_max = hm_x_min + 100.0
const hm_y_min = 0.0
const hm_y_max = 1.0
const hm_z_na = TokyoVariantAnalysis4A.lh_logL_threshold
@info "(hm_x_size, hm_y_size)", (hm_x_size, hm_y_size)
@info "(hm_x_min, hm_x_max)", (hm_x_min, hm_x_max)
@info "(hm_y_min, hm_y_max)", (hm_y_min, hm_y_max)
@info "hm_z_na", hm_z_na
heatmap_array = fill(hm_z_na, (hm_x_size, hm_y_size))
for once in 1:1
    b_array = falses(hm_x_size, hm_y_size)
    x_range = range(hm_x_min, hm_x_max; length=hm_x_size)
    for (i, row) in enumerate(eachrow(likelihoods_df))
        if i % 1000 == 0
            print(stderr, "\r$(round(Int, 100.0*i/nrow(likelihoods_df)))%")
        end
        iy_p::Union{Int,Nothing} = nothing
        for (ix, x) in enumerate(x_range)
            y = 1.0 / (1.0 + 1.0 / exp(row.α + row.β * x))
            iy = round(Int, (hm_y_size - 1.0) * (y - hm_y_min) / (hm_y_max - hm_y_min)) + 1
            if 1 ≤ iy ≤ hm_y_size
                rlogL = row.logL - likelihood_max.logL
                if !isnothing(iy_p) && (1 ≤ iy_p ≤ hm_y_size) && (iy_p != iy)
                    for iy_c in iy_p:sign(iy-iy_p):iy
                        if !b_array[ix, iy_c]
                            heatmap_array[ix, iy_c] = rlogL
                            b_array[ix, iy_c] = true
                        end
                    end
                else
                    if !b_array[ix, iy]
                        heatmap_array[ix, iy] = rlogL
                        b_array[ix, iy] = true
                    end
                end
            end
            iy_p = iy
        end
    end
    print(stderr, "\r100%\n")
end
#
function hm_f(x, y)
    ix = round(Int, (hm_x_size - 1.0) * (x - hm_x_min) / (hm_x_max - hm_x_min)) + 1
    iy = round(Int, (hm_y_size - 1.0) * (y - hm_y_min) / (hm_y_max - hm_y_min)) + 1
    if (1 ≤ ix ≤ hm_x_size) && (1 ≤ iy ≤ hm_y_size)
        return exp(heatmap_array[ix, iy])
    else
        return exp(hm_z_na)
    end
end

end #module

max_stars_repo_name: muimuiz/TokyoVariant

max_issues_repo_name: muimuiz/TokyoVariant

max_forks_repo_name: muimuiz/TokyoVariant


text:
using KnetLayers, Test, Knet, CUDA
CUDA.functional() ?  KnetLayers.settype!(KnetArray{Float64}) : KnetLayers.settype!(Array{Float64})
println("Testing with $(KnetLayers.arrtype)")


text:
using OmegaCore
using OmegaTest
using Distributions
using OmegaDistributions
using Test

function f()
  μ = 5
  xs = manynth(Normal(μ, 1), CartesianIndices((10,)))
  randsample(xs)
  # randsample(xs)
end

@testset "multivariate" begin
  f()
  μ = 5
  # xs = 1 ~ Mv(Normal(5, 1), (1000,))
  xs = manynth(Normal(μ, 1), CartesianIndices((1000,)))
  samples = randsample(xs)
  # @test isinferred(randsample, xs)
  @test mean(samples) ≈ μ atol = 0.3
end

@testset "Multivariate 2" begin
  x = 1 ~ Normal(0, 1)
  function f(id, ω)
    x(ω) + Uniform(0, 1)(id, ω)
  end
  xs = dimsnth(f, (3, 3))
  x_, xs_ = randsample((x, xs))
  @test size(xs_) == (3, 3)
  # @test isinferred(randsample, xs)
  # FIXME: not generalized to arbitrary cartesian N
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Stheno: pairwise

@testset "addition" begin
    let
        rng, N, N′, D, gpc = MersenneTwister(123456), 5, 6, 2, GPC()
        X, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        μ1, μ2, k1, k2 = ConstantMean(1.0), ConstantMean(2.0), EQ(), Linear(1.0)
        f1, f2 = GP(μ1, k1, gpc), GP(μ2, k2, gpc)
        f3 = f1 + f2
        f4 = f1 + f3
        f5 = f3 + f4

        for (n, (fp, fa, fb)) in enumerate([(f3, f1, f2), (f4, f1, f3), (f5, f3, f4)])
            Σp = cov(fa(X)) + cov(fb(X)) + xcov(fa(X), fb(X)) + xcov(fb(X), fa(X))
            ΣpXX′ = xcov(fa(X), fa(X′)) + xcov(fb(X), fb(X′)) + xcov(fa(X), fb(X′)) +
                xcov(fa(X), fb(X′))
            # @test mean(fp) == Stheno.CompositeMean(+, mean(fa), mean(fb))
            @test mean_vec(fp(X)) ≈ mean_vec(fa(X)) + mean_vec(fb(X))
            @test cov(fp(X)) ≈ Σp
            @test xcov(fp(X), fp(X′)) ≈ ΣpXX′
            @test xcov(fp(X′), fp(X)) ≈ transpose(ΣpXX′)
            @test xcov(fp(X), fa(X′)) ≈ xcov(fa(X), fa(X′)) + xcov(fb(X), fa(X′))
            @test xcov(fp(X′), fa(X)) ≈ xcov(fa(X′), fa(X)) + xcov(fb(X′), fa(X))
            @test xcov(fa(X), fp(X′)) ≈ xcov(fb(X), fa(X′)) + xcov(fa(X), fa(X′))
            @test xcov(fa(X′), fp(X)) ≈ xcov(fb(X′), fa(X)) + xcov(fa(X′), fa(X))
        end
    end

    let
        rng, N, N′, D, gpc = MersenneTwister(123456), 5, 6, 2, GPC()
        X, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        μ1, μ2, k1, k2 = ConstantMean(1.0), ConstantMean(2.0), EQ(), Linear(1.0)

        # Addition of BlockGPs.
        f1, f2 = GP(μ1, k1, gpc), GP(μ2, k2, gpc)
        fj1, fj2 = BlockGP([f1, f2]), BlockGP([f2, f1])
        @test mean_vec((fj1 + fj2)(BlockData([X, X′]))) ==
            BlockVector([mean_vec((f1 + f2)(X)), mean_vec((f2 + f1)(X′))])

        # Addition of a BlockGP and a GP.
        @test mean_vec((fj1 + f1)(BlockData([X, X′]))) ==
            BlockData([mean_vec((f1 + f1)(X)), mean_vec((f2 + f1)(X′))])

        # Addition of a GP and a BlockGP.
        @test mean_vec((f2 + fj2)(BlockData([X, X′]))) ==
            BlockData([mean_vec((f2 + f2)(X)), mean_vec((f2 + f1)(X′))])
    end

    # Check that adding a constant to a GP yields a GP with a shifted mean.
    let
        rng, N, D = MersenneTwister(123456), 5, 6, 2
        X = ColsAreObs(randn(rng, D, N))
        c, f = randn(rng), GP(5.0, EQ(), GPC())
        @test map(mean(f + c), X) == map(mean(f), X) .+ c
        @test map(mean(c + f), X) == c .+ map(mean(f), X)
        @test pairwise(kernel(f + c), X) == pairwise(kernel(f), X)
        @test pairwise(kernel(c + f), X) == pairwise(kernel(f), X)

        @test map(mean(f - c), X) == map(mean(f), X) .- c
        @test map(mean(c - f), X) == c .- map(mean(f), X)
        @test pairwise(kernel(f - c), X) == pairwise(kernel(f), X)
        @test pairwise(kernel(c - f), X) == pairwise(kernel(f), X)

        g = f(randn(rng, 10))
        @test mean_vec(g + c) == mean_vec(g) .+ c

        x = randn(rng, N + D)
        @test map(mean(f + sin), x) == map(mean(f), x) + map(sin, x)
        @test map(mean(sin + f), x) == map(sin, x) + map(mean(f), x)
        @test pairwise(kernel(f + sin), x) == pairwise(kernel(f), x)
        @test pairwise(kernel(sin + f), x) == pairwise(kernel(f), x)
    end
end


text:
push!(LOAD_PATH, "./")
# if you want to experiment with the "pushed" matrix creation method,
# switch comments on the two following lines:
using Sparse23
#using Sparse23push
using Printf
using SparseArrays,LinearAlgebra

function dotest(dim::Int,size::Int)
    # 1) matrix creation:
    t1 = time_ns()
    if dim==2
        M,order,nc=PreLapl2(size)
        sizeV=size^2
    else
        M,order,nc=PreLapl3(size)
        sizeV=size^3
    end
    t1 = time_ns()-t1
    V=Vector{Float64}([1 for i in 1:sizeV])
    W=Vector{Float64}([0 for i in 1:sizeV])

    # 2) matrix x vector product.
    t2 = time_ns()
    #A_mul_B!(W,M,V)
    mul!(W,M,V)

    t2 = time_ns()-t2

    return order,nc,t1,t2
end

println()
@printf "%-5s %1s %-20s\n" "size" ":" "size of the grid."
@printf "%-5s %1s %-20s\n" "order" ":" "order of the matrix."
@printf "%-5s %1s %-20s\n" "nc" ":" "number of ceofficients."
@printf "%-5s %1s %-20s\n" "T.b" ":" "time to build the matrix (s)."
@printf "%-5s %1s %-20s\n" "T.p" ":" "time for matrix x vector product (s)."

# computation starts here:
for dim=2:3
    println("\n"*string(dim)*"d:")
    fw=open("RunningOn"*gethostname()*"-"*string(dim),"w")
    if dim==2
        size=32
        sizeM=2048
    else
        size=16
        sizeM=256
    end
    @printf "%4s %8s %9s %10s %10s \n" "size" "order" "nc" "T.b" "T.p"
    while size<=sizeM

        order,nc,t1,t2=dotest(dim,size)#warmup

        order,nc,t1,t2=dotest(dim,size)

        s=@sprintf "%4d %8d %9d %10.2e %10.2e \n" size order nc float(t1)*10.0^(-9) float(t2)*10.0^(-9)
        print(s)
        write(fw,s)
        size*=2
    end
    close(fw)
end

max_stars_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_issues_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_forks_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo


text:
function get_probs_fv_in_clust(d::Distribution{DataSet})
  data_set = d.space
  probability_dict = Dict()
  for clust_name in keys(data_set.clusters)
    probability_dict[clust_name] = prob_clust_in_dataset(d,clust_name)
  end
  return probability_dict
end

# given an array of texts, split texts into training data and test data
function split_dataset(dataset::Array, ratio=.67)
  train_size = convert(Int64,ceil((length(dataset) * ratio)))
  train_set = Array(Any, train_size)
  dataset = shuffle!(dataset)
  i = 1
  while train_size > 0
    train_set[i] = dataset[i]
    i+=1
    train_size-=1
  end

  test_set = Array(Any, convert(Int64,((length(dataset)-i) +1)))
  j = 1
  for text in dataset[i:end]
    test_set[j] = text
    j+=1
  end

  return (train_set,test_set)
end

function separate_by_class(d::Distribution{DataSet})
  classes = Array(Any, length(d.space.clusters))
  i = 1
  for clust_name in keys(d.space.clusters)
    classes[i] = clust_name
    i+=1
  end
  return classes
end

function naive_bayes(d::Distribution{DataSet},fv::FeatureVector)

  class_probs = FeatureVector()
  for feature in keys(fv)
    for clust in keys(d.space.clusters)
      prob = cond_prob_f_given_clust(d,clust,feature)
      if isnan(prob) || prob == 0
        prob = 1e-15
      end
      class_probs[clust] += log(prob)*fv[feature]
    end
  end

  for clust in keys(d.space.clusters)
    prob = prob_clust_in_dataset(d,clust)
    class_probs[clust] += log(prob)
  end

  max_class = "unknown"
  max_value = -Inf
  num_of_same_probs = 1
  for class in keys(class_probs)
    value = class_probs[class]
    if value > max_value
      max_class = class
      max_value = value
    elseif value == max_value
      num_of_same_probs += 1
    end
  end

  if num_of_same_probs == length(d.space.clusters)
    warn("All probabilities are the same value")
  end

  return max_class
end
naive_bayes(ds::DataSet,fv::FeatureVector) = naive_bayes(Distribution(ds),fv)

function train_data(ds::DataSet,ig_list::Set)

  new_ds = DataSet()

  for clust in keys(ds.clusters)
    new_clust = Cluster()
    for fv in keys(ds[clust])
      new_fv = FeatureVector()
      for feature in keys(ds[clust][fv])
        if feature in ig_list
          new_fv[feature] = ds[clust][fv][feature]
        end
      end
      new_clust[fv] = new_fv
    end
    new_ds[clust] = new_clust
  end

  return new_ds
end

train_data(ds::DataSet,ig_list) = train_data(ds,Set(ig_list))

function percentages(correct_d::Distribution{DataSet},guesses::Dict)
  right = 0
  for guess in keys(guesses)
    if guess in keys(correct_d.space[guesses[guess]])
      right += 1
    end
  end

  correct = right/length(guesses)
  return correct
end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Reads files from
https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/
"""

using SparseArrays, LinearAlgebra

mutable struct RegressionData
  coefficient_matrix::SparseMatrixCSC{Float64,Int64}
  target_values::Vector{Float64}
end

function load_libsvm_file(file_name::String)
  open(file_name, "r") do io
    target = Array{Float64,1}()
    row_indicies = Array{Int64,1}()
    col_indicies = Array{Int64,1}()
    matrix_values = Array{Float64,1}()

    row_index = 0
    for line in eachline(io)
      row_index += 1
      split_line = split(line)
      push!(target, parse(Float64,split_line[1]))
      for i = 2:length(split_line)
        push!(row_indicies, row_index)
        matrix_coef = split(split_line[i],":")
        push!(col_indicies, parse(Int64,matrix_coef[1]))
        push!(matrix_values, parse(Float64,matrix_coef[2]))
      end
    end
    coefficient_matrix = sparse(row_indicies, col_indicies, matrix_values)
    return RegressionData(coefficient_matrix, target)
  end
end

function normalize_columns(coefficient_matrix::SparseMatrixCSC{Float64,Int64})
  m = size(coefficient_matrix,2)
  normalize_columns_by = ones(m)
  for j = 1:m
    col_vals = coefficient_matrix[:,j].nzval
    if length(col_vals) > 0
      normalize_columns_by[j] = 1.0 / norm(col_vals, 2)
    end
  end
  return coefficient_matrix * sparse(1:m,1:m,normalize_columns_by)
end

function remove_empty_columns(coefficient_matrix::SparseMatrixCSC{Float64,Int64})
  keep_cols = Array{Int64,1}()
  for j = 1:size(coefficient_matrix,2)
    if length(coefficient_matrix[:,j].nzind) > 0
      push!(keep_cols, j)
    end
  end
  return coefficient_matrix[:,keep_cols]
end

function add_intercept(coefficient_matrix::SparseMatrixCSC{Float64,Int64})
  return [sparse(ones(size(coefficient_matrix,1))) coefficient_matrix]
end

function preprocess_regression_problem(result::RegressionData)
  result.coefficient_matrix = remove_empty_columns(result.coefficient_matrix)
  result.coefficient_matrix = add_intercept(result.coefficient_matrix)
  result.coefficient_matrix = normalize_columns(result.coefficient_matrix)
  return result
end


text:
#using Pkg; Pkg.add("LinearAlgebra"); Pkg.add("Markdown")
# using Documenter
using LinearAlgebra
using Markdown                             # Pour que les docstrings en début des fonctions ne posent
                                           # pas de soucis. Ces docstrings sont utiles pour générer
                                           # la documentation sous GitHub
include("Algorithme_De_Newton.jl")

# Affichage les sorties de l'algorithme des Régions de confiance
function my_afficher_resultats(algo,nom_fct,point_init,xmin,fxmin,flag,sol_exacte,nbiters)
	println("-------------------------------------------------------------------------")
	printstyled("Résultats de : ",algo, " appliqué à ",nom_fct, " au point initial ", point_init, ":\n",bold=true,color=:blue)
	println("  * xsol = ",xmin)
	println("  * f(xsol) = ",fxmin)
	println("  * nb_iters = ",nbiters)
	println("  * flag = ",flag)
	println("  * sol_exacte : ", sol_exacte)
end

# Fonction f0
# -----------
f0(x) =  sin(x)
# la gradient de la fonction f0
grad_f0(x) = cos(x)
# la hessienne de la fonction f0
hess_f0(x) = -sin(x)
sol_exacte = -pi/2
options = []

x0 = sol_exacte
xmin,f_min,flag,nb_iters = Algorithme_De_Newton(f0,grad_f0,hess_f0,x0,options)
my_afficher_resultats("Newton","f0",x0,xmin,f_min,flag,sol_exacte,nb_iters)
x0 = -pi/2+0.5
xmin,f_min,flag,nb_iters = Algorithme_De_Newton(f0,grad_f0,hess_f0,x0,options)
my_afficher_resultats("Newton","f0",x0,xmin,f_min,flag,sol_exacte,nb_iters)
x0 = pi/2
xmin,f_min,flag,nb_iters = Algorithme_De_Newton(f0,grad_f0,hess_f0,x0,options)
my_afficher_resultats("Newton","f0",x0,xmin,f_min,flag,sol_exacte,nb_iters)

max_stars_repo_name: faicaltoubali/ENSEEIHT

max_issues_repo_name: faicaltoubali/ENSEEIHT

max_forks_repo_name: faicaltoubali/ENSEEIHT


text:
"""
    cleanWrite()

Clears all timestamp directories in the current folder
"""
function cleanWrite()
    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
    for i =1:length(dirresults)
        rm("$(dirresults[i])", force=true, recursive=true)
    end
    rm("*~", force=true)
end

"""
    simpleInterp(x1, x2, y1, y2, x)

Performs a linear interpolation between points (x1,y1) and (x2,y2) to
find the value of `y` at `x`.

"""
function simpleInterp(x1 ::Float64, x2 :: Float64, y1 :: Float64, y2 :: Float64, x::Float64)
    y = y1 + (y2 - y1)*(x - x1)/(x2 - x1)
    return y
end

"""
    find_tstep([kin])

kin is a definition of either Sin / Cos / Ramp kinematics

A suitable time-step is calculated based on the values of amplitude and frequency

For Sin / Cos kinematics, time step for amp*k = 0.2 is 0.015

For Eldredge ramp kinematics, time step for K = 0.2 is 0.015

For Constant kinematics, time step is 0.015

"""
function find_tstep(kin:: CosDef)
    dtstar = 15.
    dt_tmp = 0.015*0.2/(kin.k*kin.amp)
    dtstar = minimum([dt_tmp dtstar])
    return dtstar
end

function find_tstep(kin:: SinDef)
    dtstar = 15.
    dt_tmp = 0.015*0.2/(kin.k*kin.amp)
    dtstar = minimum([dt_tmp dtstar])
    return dtstar
end

function find_tstep(kin:: EldUpDef)
    dtstar = 1.
    dtstar = minimum([0.015*0.2/kin.K 0.015])
    return dtstar
end

function find_tstep(kin:: EldUptstartDef)
    dtstar = 1.
    dtstar = minimum([0.015*0.2/kin.K 0.015])
    return dtstar
end

function find_tstep(kin:: EldRampReturnDef)
    dtstar = 1.
    dtstar = minimum([0.015*0.2/kin.K 0.015])
    return dtstar
end

function find_tstep(kin :: BendingDef)
    dtstar = 15
    amp = evaluate(kin.spl, kin.spl.t[end])
    dt_tmp = 0.015*0.2/(kin.k*amp)
    dtstar = minimum([dt_tmp dtstar])
    return dtstar
end

function find_tstep(kin :: ConstDef)
    dtstar = 0.015
    return dtstar
end

# Numerical integration method: Trapezoidal
function simpleTrapz(y::Vector{T}, x::Vector{T}) where {T<:Real}
    local len = length(y)
    if len != length(x)
        error("Vectors must be of same length")
    end
    r = 0.0
    for i in 2:len
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    r/2.0
end

# Aerofoil camber calculation from coordinate file
function camber_calc(x::Vector,airfoil::String)
    #Determine camber and camber slope on airfoil from airfoil input file

    ndiv = length(x);
    c = x[ndiv];

    cam = zeros(ndiv)
    cam_slope = zeros(ndiv)
    in_air = DelimitedFiles.readdlm(airfoil, Float64);
    xcoord = in_air[:,1];
    ycoord = in_air[:,2];
    ncoord = length(xcoord);
    xcoord_sum = zeros(ncoord);

    for i = 1:ncoord-1
        xcoord_sum[i+1] = xcoord_sum[i] + abs(xcoord[i+1]-xcoord[i]);
    end
    y_spl = Spline1D(xcoord_sum,ycoord);

    y_ans = zeros(2*ndiv);

    for i=1:ndiv
        y_ans[i] = evaluate(y_spl,x[i]/c);
    end
    for i=ndiv+1:2*ndiv
        y_ans[i] = evaluate(y_spl,(x[ndiv]/c) + (x[i-ndiv]/c));
    end
    cam[1:ndiv] = [(y_ans[i] + y_ans[(2*ndiv) + 1 - i])*c/2 for i = ndiv:-1:1];
    cam[1] = 0;
    cam_spl = Spline1D(x,cam);
    cam_slope[1:ndiv] = derivative(cam_spl,x);
    return cam, cam_slope
end

function camber_thick_calc(x::Vector,coord_file::String)
    #Determine camber and camber slope on airfoil from airfoil input file

    ndiv = length(x);
    c = x[ndiv];

    cam = zeros(ndiv)
    cam_slope = zeros(ndiv)
    thick = zeros(ndiv)
    thick_slope = zeros(ndiv)

    if coord_file[1:6] == "NACA00"
        m = parse(Int, coord_file[5])/100.
        p = parse(Int, coord_file[6])/10.
        th = parse(Int, coord_file[7:8])/100.

        b1 = 0.2969; b2 = -0.1260; b3 = -0.3516; b4 = 0.2843; b5 = -0.1015
        for i = 2:ndiv
            thick[i] = 5*th*(b1*sqrt(x[i]) + b2*x[i] + b3*x[i]^2 + b4*x[i]^3 + b5*x[i]^4)
            thick_slope[i] = 5*th*(b1/(2*sqrt(x[i])) + b2 + 2*b3*x[i] + 3*b4*x[i]^2 + 4*b5*x[i]^3)
        end
        thick[1] = 5*th*(b1*sqrt(x[1]) + b2*x[1] + b3*x[1]^2 + b4*x[1]^3 + b5*x[1]^4)
        thick_slope[1] = 2*thick_slope[2] - thick_slope[3]
        rho = 1.1019*th*th*c
        cam[1:ndiv] .= 0.
        cam_slope[1:ndiv] .= 0.

    elseif coord_file[1:8] == "Cylinder"
        for i = 1:ndiv
            theta = acos(1. - 2*x[i]/c)
            thick[i] = 0.5*c*sin(theta)
        end
        thick_spl = Spline1D(x,thick)
        thick_slope[1:ndiv] = derivative(thick_spl,x)
        cam[1:ndiv] .= 0.
        cam_slope[1:ndiv] .= 0.
        rho = 0.5*c
    elseif coord_file[1:9] == "FlatPlate"
        th = parse(Int, coord_file[10:13])/10000.
        r = th*c/2
        for i = 2:ndiv-1
            if x[i] <= r
                thick[i] = sqrt(r^2 - (x[i] - r)^2)
                thick_slope[i] = -(x[i] - r)/(sqrt(r^2 - (x[i] - r)^2))
            elseif  x[i] >= c-r
                thick[i] = sqrt(r^2 - (x[i] - c + r)^2)
                thick_slope[i] = -(x[i] - c + r)/sqrt(r^2 - (x[i] - c + r)^2)
            else
                thick[i] = r
            end
        end
        thick[1] = sqrt(r^2 - (x[1] - r)^2)
        thick_slope[1] = 2*thick_slope[2] - thick_slope[3]
        thick[ndiv] = sqrt(r^2 - (x[ndiv] - c + r)^2)
        thick_slope[ndiv] = 2*thick_slope[ndiv-1] - thick_slope[ndiv-2]

        rho = r
        for i = 1:ndiv
            cam[i] = 0.
            cam_slope[i] = 0.
        end
    else
        coord = readdlm(coord_file)
        ncoord = length(coord[:,1])
        if (0. in coord[:,1]) == false
            error("Airfoil file must contain leading edge coordinate (0,0)")
        else
            nle = find(x->x==0, coord[:,1])[1]
            if coord[nle,2] != 0.
                error("Airfoil leading edge must be at (0,0)")
            end
        end

        zu_spl = Spline1D(reverse(coord[1:nle,1]), reverse(coord[1:nle,2]),k=1)
        zl_spl = Spline1D(coord[nle:ncoord,1], coord[nle:ncoord,2],k=1)

        zu = Array{Float64}(ndiv)
        zl = Array{Float64}(ndiv)

        for i=1:ndiv
            zu[i] = evaluate(zu_spl,x[i]/c)
            zl[i] = evaluate(zl_spl,x[i]/c)
        end

        cam[1:ndiv] = [(zu[i] + zl[i])*c/2 for i = 1:ndiv]
        thick[1:ndiv] = [(zu[i] - zl[i])*c/2 for i = 1:ndiv]
        cam_spl = Spline1D(x,cam)
        thick_spl = Spline1D(x,thick)
        cam_slope[1:ndiv] = derivative(cam_spl,x)
        thick_slope[1:ndiv] = derivative(thick_spl,x)
        rho = readdlm("rho")[1]
    end
    return thick, thick_slope,rho, cam, cam_slope
end

max_stars_repo_name: hjabird/UNSflow

max_issues_repo_name: hjabird/UNSflow

max_forks_repo_name: hjabird/UNSflow


text:
module IsopycnalSurfaces

using Dierckx, Interpolations

export sigma0column, sigma1column, sigma2column,
 vars2sigma1,  var2sigmacolumn, sigma1grid,
 mixinversions!, dedup!

"""
    function vars2sigma1(vars,p,sig1grid,γ,spline_order)
    map variables from regular 3D grid onto sigma1 surfaces
# Arguments
- `vars::Dict{String,Array{T,3}}}`: dict of 3d arrays
- `p::Array{T,1}` : vertical profile of standard pressures
- `sig1grid`: σ₁ surface values
- `γ`: grid description needed for preallocation
- `splorder`: 1-5, order of spline
# Output
- `varsσ::Dict{String,Array{T,3}`: dict of 3d arrays of variables on sigma1 surfaces
"""
function vars2sigma1(vars::Dict{String,Array{T,3}},pressure::Vector{T},σ₁grid::Vector{T},splorder::Integer,linearinterp=false) where T<:AbstractFloat

    # is there a problem if pressure is not the same type as the input vars?
    # could introduce two parametric types to function definition above

    # θ and S must exist

    # a list of potential temperature names that are recognized
    θkeys = ("THETA","θ","theta","T")
    Skeys = ("SALT","S","Sp")

    θname = string()

    nx = Integer(); ny = Integer(); nz = Integer()

    for k in θkeys
        if haskey(vars,k)
            nx,ny,nz = size(vars[k])
            θname = k
        end
    end
    θname == "" && error("Potential temperature missing")

    Sname = string()
    for k in Skeys
        if haskey(vars,k)
            nx,ny,nz = size(vars[k])
            Sname = k
        end
    end
    Sname == "" && error("Salinity missing")

    # loop over faces
    nσ = length(σ₁grid)

    # vcol = Dict with profile/column data
    # pre-allocate each key in vars
    vcol = Dict{String,Array{T,1}}() # vars in a column
    varsσ = Dict{String,Array{T,3}}()

    for (key, value) in vars
        vcol[key] = fill(convert(T,NaN),nz)
        varsσ[key] = fill(convert(T,NaN),(nx,ny,nσ))
    end

    # allocate standard pressure by hand.
    varsσ["p"] = fill(convert(T,NaN),(nx,ny,nσ))

    for xx = 1:nx
        for yy = 1:ny
            for (vcolname, vcolval) in vars
                # vcol = Dict with profile/column data
                vcol[vcolname] = vcolval[xx,yy,:]
            end

            # also need to filter dry values and to change zz
            # Consider using `isdry` function and dryval in future.
            nw = count(notnanorzero,vcol[θname]) # number of wet points in column
            nwS = count(notnanorzero,vcol[Sname])
            if nw != nwS
                error("T,S zeroes inconsistent")
            end

            if nw > 3
                # incurs error if splorder > number of points in column
                # if nw > splorder #need >=n+1 points to do order-n interpolation
                σ₁=sigma1column(vcol[θname][1:nw],vcol[Sname][1:nw],pressure[1:nw])

                for (vckey,vcval) in vcol
                    varσ = var2sigmacolumn(σ₁,vcval[1:nw],σ₁grid,splorder,linearinterp)
                    [varsσ[vckey][xx,yy,ss] = convert(T,varσ[ss]) for ss = 1:nσ]
                end

                # do standard pressure by hand.
                pσ = var2sigmacolumn(σ₁,pressure[1:nw],σ₁grid,splorder,linearinterp)
                [varsσ["p"][xx,yy,ss] = convert(T,pσ[ss]) for ss = 1:nσ]

            end
        end
    end
    return varsσ
end

"""
    function sigmacolumn(θ,S,p,p0)
    σ for a water column
# Arguments
- `θz::Vector{T}`: potential temperature
- `Sz::Vector{T}`: practical salinity
- `pz::Vector{T}`: vertical profile of standard pressures
- `p0::Integer`: reference pressure
# Output
- `σ::Vector{T}`:  sigma for wet points in column
"""
function sigmacolumn(θz::Vector{T},Sz::Vector{T},pz::Vector{T2},p0::Integer)::Vector{T} where T<:AbstractFloat where T2<:AbstractFloat
    nz = length(θz)
    σ = similar(θz)

    σa,σb,σc = SeaWaterDensity(θz,Sz,pz,p0)
    [σ[zz] = convert(T,σc[zz]) .- 1000.0 for zz = 1:nz]
    return σ
end

"""
    function sigma0column(θ,S,p)
    σ₀ for a water column
    Untested for a mix of float values

# Arguments
- `θz::Vector{T}`: potential temperature
- `Sz::Vector{T}`: practical salinity
- `pz::Vector{T}`: vertical profile of standard pressures
# Output
- `σ₀`:  sigma-0 for wet points in column
"""
sigma0column(θz,Sz,pz) = sigmacolumn(θz,Sz,pz,0)

"""
    function sigma1column(θ,S,p)
    σ₁ for a water column
    Untested for a mix of float values

# Arguments
- `θz::Vector{T}`: potential temperature
- `Sz::Vector{T}`: practical salinity
- `pz::Vector{T}`: vertical profile of standard pressures
# Output
- `σ₁`:  sigma-1 for wet points in column
"""
sigma1column(θz,Sz,pz) = sigmacolumn(θz,Sz,pz,1000)

"""
    function sigma2column(θ,S,p)
    σ₂ for a water column
    Untested for a mix of float values

# Arguments
- `θz::Vector{T}`: potential temperature
- `Sz::Vector{T}`: practical salinity
- `pz::Vector{T}`: vertical profile of standard pressures
# Output
- `σ₂`:  sigma-2 for wet points in column
"""
sigma2column(θz,Sz,pz) = sigmacolumn(θz,Sz,pz,2000)

"""
   function notnanorzero

     true is argument is not a NaN nor zero
"""
notnanorzero(z) = !iszero(z) && !isnan(z)

""" function sigma1grid()
    Standard (from Susan Wijffels, WHOI) choice of sigma1 surfaces
# Arguments
- `z`: value
# Output
- `σ₁grid`: list (vector) of σ₁ values
"""
function sigma1grid()
    σ₁grida = 24:0.05:31
    σ₁gridb = 31.02:0.02:33
    σ₁grid = vcat(σ₁grida,σ₁gridb)
    σ₁grid = σ₁grid[1:3:end]
    return σ₁grid
end

"""
    function var2sigmacolumn(σ,v,σ₁grid,splorder)
    map θ,S, p onto σ₁ surfaces for a water column
# Arguments
- `σ`::Array{Float,1}}`: sigma values of input variable
- `v::Array{Float,1}}`: variable of interest
- `sig1`: σ surface values
- `splorder`: 1-5, order of spline
- `linearinterp`: optional argument, true to force linear interpolation
# Output
- `θonσ`: variable on sig1 sigma surfaces
"""
function var2sigmacolumn(σorig::Vector{T},v,σgrid,splorder::Integer,linearinterp=false) where T<:AbstractFloat where T2<:AbstractFloat
    # choose a univariate spline with s = magic number
    #θspl = Spline1D(σ₁,θz;k=splorder,s=length(σ₁))

    σ = copy(σorig) # make sure sigma-1 doesn't mutate and pass back

    nσout = length(σgrid)
    θonσ = fill(convert(T,NaN),nσout)

    # some velocity variables may include NaN's
    # add this kludge to make NaN velocity equal to zero
    # a useful approximation to no-flow boundary condition?
    replace!(v,convert(T,NaN)=>zero(T))

    # eliminate homogeneities
    dedup!(σ,v)

    # mix inversions
    mixinversions!(σ,v)

    nσin = length(σ)

    # 1) no inversions or homogeneity (this constraint relaxed now because too many profiles thrown out)
    # 2) range of sig1, 3) no extrapolation
    # if sum(diff(σ₁).<0)==0 && count(minimum(σ₁).<=σ₁grid.<=maximum(σ₁)) > 0
    # added nσin > 1 to fix error
    if nσin > 1 && count(minimum(σ) .<= σgrid .<= maximum(σ)) > 0

        # eliminate any extrapolation
        sgood = findall(minimum(σ).<=σgrid.<= maximum(σ))

        if nσin > splorder
            θspl = Spline1D(σ,v;k=splorder)
            for ss in sgood
                θonσ[ss] = θspl(σgrid[ss])
            end

            # check for spline instability

            # this version didn't work if length(sgood)==1
            #if maximum(θonσ[sgood]) - minimum(θonσ[sgood]) > 1.00 * (maximum(v) - minimum(v))

            # give some leeway with "0.1"
            if (maximum(θonσ[sgood]) > maximum(v) + 0.1*(maximum(v)-minimum(v)) ||
                minimum(θonσ[sgood]) < minimum(v) - 0.1*(maximum(v)-minimum(v)))
                linearinterp = true
                #println("unstable spline")
            end
        else # spline interp
            linearinterp = true
        end

        if linearinterp
            #println("doing linear interp")
#            println(size(σ),size(v))
            interp_linear = LinearInterpolation(σ, v)
            for ss in sgood
                θonσ[ss] = interp_linear(σgrid[ss])
            end
        end # linearinterp

    else
        #println("not doing any calcs")
    end # any good points?

    return θonσ
end

"""
function mixinversions!(a,b)

For values of `a` that are not increasing,
an inversion is defined.
Average these values of a until values of `a` are sorted.
Do the same averaging on the accompanying vector `b`.

# Arguments
- `a::Vector{T}`: a density variable
- `b::Vector{T}`: an accompanying tracer variable

Arguments are mutated by the function.
"""
function mixinversions!(a,b)
    while sum(diff(a).<=0) > 0
        length(a) == 1 ? da = 1. : da = diff(a)
        ii = findfirst(isnotpositive,da)
        a[ii] = (a[ii] + a[ii+1])/2
        b[ii] = (b[ii] + b[ii+1])/2
        deleteat!(a,ii+1)
        deleteat!(b,ii+1)
    end
end

"""
function dedup!(a,b)

Remove values of `a` that are duplicates.
Remove values of `b` that have the same location as the duplicates in `a`.
The length of `a` and `b` should be identical before and after invoking this function.

# Arguments
- `a::Vector{T}`: a density variable
- `b::Vector{T}`: an accompanying tracer variable

Arguments are mutated by the function.
"""
function dedup!(a,b)
    length(a) == 1 ? da = 1. : da = diff(a)
    while count(iszero,da) > 0
        dedupfirst!(a,b)
        length(a) ==1 ? da = 1. : da = diff(a)
    end
end

"""
function dedup!(a,b)

Remove the first duplicate of `a`.
Remove the value of `b` that is located at the same entry as the first duplicate in `a`.
The length of `a` and `b` should be identical before and after invoking this function.

# Arguments
- `a::Vector{T}`: a density variable
- `b::Vector{T}`: an accompanying tracer variable

Arguments are mutated by the function.
"""
function dedupfirst!(a,b)
    counter = 1
    da = diff(a) # requires length of 2 or more
    ii = findfirst(iszero,da)
    while iszero(da[ii])
        b[ii+1] += b[ii]
        deleteat!(a,ii)
        deleteat!(b,ii)
        deleteat!(da,ii)
        counter += 1
        (length(da) < ii) ? break : nothing
    end
    b[ii] /= counter
end

isnotpositive(x) = (abs(x) == -x)

"""
SeaWaterDensity(Θ,Σ,Π,Π0) from MITgcmTools.jl/PhysicalOceanography.jl, From Gael Forget

Compute potential density (ρP), in situ density (ρI), and density
referenced to PREF (Π0 in decibars) from potential temperature (Θ in °C),
salinity (Σ in psu) and pressure (Π in decibars) according to the
UNESCO / Jackett & McDougall 1994 equation of state.

Credits: code based on a Matlab implementation by B. Ferron

Reference: https://www.jodc.go.jp/info/ioc_doc/UNESCO_tech/059832eb.pdf

Check value: ρI = `1041.83267kg/m^3` for Θ=`3°Celcius`, Σ=`35psu`, Π=`3000dbar`
```
(ρP,ρI,ρR) = SeaWaterDensity(3.,35.5,3000.)
isapprox(ρI,1041.83267, rtol=1e-6)
```
"""
function SeaWaterDensity(Θ,Σ,Π,Π0=missing)

   #square root salinity
   sqrtΣ= sqrt.(Σ)
   #compute density pure water at atm pressure
   ZR1= ((((6.536332E-9*Θ .-1.120083E-6).*Θ .+1.001685E-4).*Θ
   .-9.095290E-3).*Θ .+6.793952E-2).*Θ .+999.842594
   #seawater density atm pressure
   ZR2= (((5.3875E-9*Θ .-8.2467E-7).*Θ .+7.6438E-5).*Θ
   .-4.0899E-3).*Θ .+0.824493
   ZR3= (-1.6546E-6*Θ .+1.0227E-4).*Θ .-5.72466E-3
   ZR4= 4.8314E-4

   #potential density (referenced to the surface)
   ρP= (ZR4*Σ + ZR3.*sqrtΣ + ZR2).*Σ + ZR1

   #add the compression terms
   ZE = (-3.508914E-8*Θ .-1.248266E-8).*Θ .-2.595994E-6
   ZBW= ( 1.296821E-6*Θ .-5.782165E-9).*Θ .+1.045941E-4
   ZB = ZBW + ZE .* Σ

   ZD = -2.042967E-2
   ZC = (-7.267926E-5*Θ .+2.598241E-3).*Θ .+0.1571896
   ZAW= ((5.939910E-6*Θ .+2.512549E-3).*Θ .-0.1028859).*Θ .-4.721788
   ZA = ( ZD*sqrtΣ + ZC).*Σ + ZAW

   ZB1= (-0.1909078*Θ .+7.390729).*Θ .-55.87545
   ZA1= ((2.326469E-3*Θ .+1.553190).*Θ .-65.00517).*Θ .+1044.077
   ZKW= (((-1.361629E-4*Θ .-1.852732E-2).*Θ .-30.41638).*Θ
   .+2098.925).*Θ .+190925.6
   ZK0= (ZB1.*sqrtΣ + ZA1).*Σ + ZKW

   #in situ density
   ρI = ρP ./ (1.0 .-Π./(ZK0-Π.*(ZA-Π.*ZB)))

   #density referenced to level Π0
   if !ismissing(Π0)
      ρR = ρP ./ (1.0 .-Π0./(ZK0-Π0.*(ZA-Π0.*ZB)))
   else
      ρR = ρP
   end

   return ρP,ρI,ρR
end

end


text:
#=
Code for the preprint "Neural ODE for Reinforcement Learning and Nonlinear Optimal Control: Cartpole Problem Revisited"
https://medium.com/@paulshen_62280/neural-ode-for-reinforcement-learning-and-nonlinear-optimal-control-cartpole-problem-revisited-5408018b8d71

Hosted on Google Colab notebook https://colab.research.google.com/drive/1p2cg_7SNG-YMhlV2mAI_53Xvbhyka4FK?usp=sharing

Paul Shen
pxshen@alumni.stanford.edu
=#

using DiffEqFlux,
    DifferentialEquations,
    Flux,
    Plots,
    Interpolations,
    DataFrames,
    CSV,
    JSON,
    Plots,
    Dates,
    Dierckx,
    FiniteDifferences,
    Optim

# physical params
m = 1 # pole mass kg
M = 2 # cart mass kg
L = 1 # pole length m
g = 9.8 # acceleration constant m/s^2

# map angle to [-pi, pi)
modpi(theta) = mod2pi(theta + pi) - pi

#=
system dynamics derivative

du: du/dt, state vector derivative updated inplace
u: state vector (x, dx, theta, dtheta)
p: parameter function, here lateral force exerted by cart as a fn of time
t: time
=#
function cartpole(du, u, p, t)
  # position (cart), velocity, pole angle, angular velocity
  x, dx, theta, dtheta = u
  force = p(t)

  du[1] = dx
  du[2] =
    (force + m * sin(theta) * (L * dtheta^2 - g * cos(theta))) /
    (M + m * sin(theta)^2)
  du[3] = dtheta
  du[4] =
    (
      -force * cos(theta) - m * L * dtheta^2 * sin(theta) * cos(theta) +
      (M + m) * g * sin(theta)
    ) / (L * (M + m * sin(theta)^2))
end

# neural network controller, here a simple MLP
# inputs: cos(theta), sin(theta) [to avoid discontinuity in theta mod 2pi], theta_dot
# output: cart force
controller = FastChain((x, p) -> x, FastDense(3, 8, tanh), FastDense(8, 1))

# initial neural network weights
pinit = initial_params(controller)

#=
system dynamics derivative with the controller included
=#
function cartpole_controlled(du, u, p, t)
  # controller force response
  force = controller([cos(u[3]), sin(u[3]), u[4]], p)[1]
  du[5] = force

  # plug force into system dynamics
  cartpole(du, u[1:4], t -> force, t)
end

# initial condition
u0 = [0; 0; pi; 0; 0]
tspan = (0.0, 1.0)
N = 50
tsteps = range(tspan[1], length = N, tspan[2])
dt = (tspan[2] - tspan[1]) / N

# set up ODE problem
prob = ODEProblem(cartpole_controlled, u0, tspan, pinit)

# wrangles output from ODE solver
function format(pred)
  x = pred[1, :]
  dx = pred[2, :]

  theta = modpi.(pred[3, :])
  dtheta = pred[4, :]

  # take derivative of impulse to get force
  impulse = pred[5, :]
  tmp = (impulse .- circshift(impulse, 1)) / dt
  force = [tmp[2], tmp[2:end]...]

  return x, dx, theta, dtheta, force
end

# solves ODE
function predict_neuralode(p)
  tmp_prob = remake(prob, p = p)
  solve(tmp_prob, Tsit5(), saveat = tsteps)
end

# loss to minimize as a function of neural network parameters p
function loss_neuralode(p)
  pred = predict_neuralode(p)
  x, dx, theta, dtheta, force = format(pred)
  loss = sum(theta .^ 2) / N + 4theta[end]^2 + dx[end]^2

  return loss, pred
end

i = 0 # training epoch counter
data = 0 # time series of state vector and control signal
# callback function after each training epoch
callback = function (p, l, pred; doplot = true)
  global i += 1

  global data = format(pred)
  x, dx, theta, dtheta, force = data

  # ouput every few epochs
  if i % 50 == 0
    println(l)
    display(plot(tsteps, theta))
    display(plot(tsteps, x))
    display(plot(tsteps, force))
  end

  return false

end

result = DiffEqFlux.sciml_train(
  loss_neuralode,
  pinit,
  ADAM(0.05),
  cb = callback,
  maxiters = 1500,
)

p = result.minimizer

# save model and data
open(io -> write(io, json(p)), "model.json", "w")
open(io -> write(io, json(data)), "data.json", "w")

gr()
x, dx, theta, dtheta, force = data
anim = Animation()

plt=plot(tsteps,[modpi.(theta.+.01),x,force],title=["Angle" "Position" "Force"],layout=(3,1))
display(plt)
savefig(plt,"cartpole_data.png")

for (x, theta) in zip(x, theta)

    cart = [x - 1 x + 1; 0 0]
    pole = [x x + 10*sin(theta); 0 10*cos(theta)]
    plt = plot(
        cart[1, :],
        cart[2, :],
        xlim = (-10, 10),
        ylim = (-10, 10),
        title = "Cartpole",
        linewidth = 3,
    )
    plot!(plt, pole[1, :], pole[2, :], linewidth = 6)

    frame(anim)
end

gif(anim, "cartpole_animation.gif", fps = 10)

max_stars_repo_name: matsuren/neural-ode-cartpole

max_issues_repo_name: matsuren/neural-ode-cartpole

max_forks_repo_name: matsuren/neural-ode-cartpole


text:
int_rules_4_1_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.5*trig^m*(a*cos+b*sin)^n =#
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~a * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n) / (~b * ~d * ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-1 / ~d) * Subst(Antiderivative((((~a) ^ 2 + (~b) ^ 2) - (~x) ^ 2) ^ ((~n - 1) / 2), ~x), ~x, ~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && IGtQ((~n - 1) / 2, 0))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-((~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x))) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * ~n) + (((~n - 1) * ((~a) ^ 2 + (~b) ^ 2)) / ~n) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(IntegerQ((~n - 1) / 2)) && GtQ(~n, 1)))
    @apply_utils Antiderivative(1 / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (-1 / ~d) * Subst(Antiderivative(1 / (((~a) ^ 2 + (~b) ^ 2) - (~x) ^ 2), ~x), ~x, ~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ 2, ~x) => sin(~c + ~d * ~x) / (~a * ~d * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1)) / (~d * (~n + 1) * ((~a) ^ 2 + (~b) ^ 2)) + ((~n + 2) / ((~n + 1) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((~a) ^ 2 + (~b) ^ 2) ^ (~n / 2) * Antiderivative(cos((~c + ~d * ~x) - atan(~a, ~b)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (Not(GeQ(~n, 1) || LeQ(~n, -1)) && GtQ((~a) ^ 2 + (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n / ((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) / sqrt((~a) ^ 2 + (~b) ^ 2)) ^ ~n) * Antiderivative(cos((~c + ~d * ~x) - atan(~a, ~b)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (Not(GeQ(~n, 1) || LeQ(~n, -1)) && Not(GtQ((~a) ^ 2 + (~b) ^ 2, 0) || EqQ((~a) ^ 2 + (~b) ^ 2, 0)))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~a) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1) * sin(~c + ~d * ~x) ^ (~n - 1)) + 2 * ~b * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1) / sin(~c + ~d * ~x) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~m + ~n, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~n, 1)))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1) * cos(~c + ~d * ~x) ^ (~n - 1)) + 2 * ~a * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1) / cos(~c + ~d * ~x) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~m + ~n, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~n, 1)))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~a * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n) / (2 * ~b * ~d * ~n * sin(~c + ~d * ~x) ^ ~n) + (1 / (2 * ~b)) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1) / sin(~c + ~d * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~m + ~n, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, 0)))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~b) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n) / (2 * ~a * ~d * ~n * cos(~c + ~d * ~x) ^ ~n) + (1 / (2 * ~a)) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1) / cos(~c + ~d * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~m + ~n, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, 0)))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((~a * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n) / (2 * ~b * ~d * ~n * sin(~c + ~d * ~x) ^ ~n)) * HypergeometricFunctions._₂F₁(1, ~n, ~n + 1, (~b + ~a * cot(~c + ~d * ~x)) / (2 * ~b)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~m + ~n, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && Not(IntegerQ(~n))))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((-(~b) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n) / (2 * ~a * ~d * ~n * cos(~c + ~d * ~x) ^ ~n)) * HypergeometricFunctions._₂F₁(1, ~n, ~n + 1, (~a + ~b * tan(~c + ~d * ~x)) / (2 * ~a)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~m + ~n, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && Not(IntegerQ(~n))))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => Antiderivative((~b + ~a * cot(~c + ~d * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~m + ~n, 0) && (IntegerQ(~n) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => Antiderivative((~a + ~b * tan(~c + ~d * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~m + ~n, 0) && (IntegerQ(~n) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (1 / ~d) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * ~x) ^ ~n) / (1 + (~x) ^ 2) ^ ((~m + ~n + 2) / 2), ~x), ~x, tan(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~n) && (IntegerQ((~m + ~n) / 2) && (NeQ(~n, -1) && Not(GtQ(~n, 0) && GtQ(~m, 1)))))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-1 / ~d) * Subst(Antiderivative(((~x) ^ ~m * (~b + ~a * ~x) ^ ~n) / (1 + (~x) ^ 2) ^ ((~m + ~n + 2) / 2), ~x), ~x, cot(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~n) && (IntegerQ((~m + ~n) / 2) && (NeQ(~n, -1) && Not(GtQ(~n, 0) && GtQ(~m, 1)))))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandTrig(sin(~c + ~d * ~x) ^ ~m * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~m) && IGtQ(~n, 0))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandTrig(cos(~c + ~d * ~x) ^ ~m * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~m) && IGtQ(~n, 0))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~a) ^ ~n * (~b) ^ ~n * Antiderivative(sin(~c + ~d * ~x) ^ ~m * (~b * cos(~c + ~d * ~x) + ~a * sin(~c + ~d * ~x)) ^ -(~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~a) ^ ~n * (~b) ^ ~n * Antiderivative(cos(~c + ~d * ~x) ^ ~m * (~b * cos(~c + ~d * ~x) + ~a * sin(~c + ~d * ~x)) ^ -(~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n / sin(~(c') + ~(d') * ~x), ~x) => (~a * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1)) + ~b * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1), ~x) + (~a) ^ 2 * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2) / sin(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n / cos(~(c') + ~(d') * ~x), ~x) => (-(~b) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1)) + ~a * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1), ~x) + (~b) ^ 2 * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2) / cos(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => -(((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(sin(~c + ~d * ~x) ^ (~m + 2) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2), ~x) + 2 * ~b * Antiderivative(sin(~c + ~d * ~x) ^ (~m + 1) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1), ~x) + (~a) ^ 2 * Antiderivative(sin(~c + ~d * ~x) ^ ~m * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (GtQ(~n, 1) && LtQ(~m, -1)))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => -(((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ (~m + 2) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2), ~x) + 2 * ~a * Antiderivative(cos(~c + ~d * ~x) ^ (~m + 1) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 1), ~x) + (~b) ^ 2 * Antiderivative(cos(~c + ~d * ~x) ^ ~m * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (GtQ(~n, 1) && LtQ(~m, -1)))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (~b * ~x) / ((~a) ^ 2 + (~b) ^ 2) - (~a / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (~a * ~x) / ((~a) ^ 2 + (~b) ^ 2) + (~b / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~m / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (-(~a) * sin(~c + ~d * ~x) ^ (~m - 1)) / (~d * ((~a) ^ 2 + (~b) ^ 2) * (~m - 1)) + (~b / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(sin(~c + ~d * ~x) ^ (~m - 1), ~x) + ((~a) ^ 2 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(sin(~c + ~d * ~x) ^ (~m - 2) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~m, 1))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~m / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (~b * cos(~c + ~d * ~x) ^ (~m - 1)) / (~d * ((~a) ^ 2 + (~b) ^ 2) * (~m - 1)) + (~a / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ (~m - 1), ~x) + ((~b) ^ 2 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ (~m - 2) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~m, 1))
    @apply_utils Antiderivative(1 / (sin(~(c') + ~(d') * ~x) * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x))), ~x) => (1 / ~a) * Antiderivative(cot(~c + ~d * ~x), ~x) - (1 / ~a) * Antiderivative((~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / (cos(~(c') + ~(d') * ~x) * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x))), ~x) => (1 / ~b) * Antiderivative(tan(~c + ~d * ~x), ~x) + (1 / ~b) * Antiderivative((~b * cos(~c + ~d * ~x) - ~a * sin(~c + ~d * ~x)) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~m / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (sin(~c + ~d * ~x) ^ (~m + 1) / (~a * ~d * (~m + 1)) - (~b / (~a) ^ 2) * Antiderivative(sin(~c + ~d * ~x) ^ (~m + 1), ~x)) + (((~a) ^ 2 + (~b) ^ 2) / (~a) ^ 2) * Antiderivative(sin(~c + ~d * ~x) ^ (~m + 2) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~m / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (-(cos(~c + ~d * ~x) ^ (~m + 1)) / (~b * ~d * (~m + 1)) - (~a / (~b) ^ 2) * Antiderivative(cos(~c + ~d * ~x) ^ (~m + 1), ~x)) + (((~a) ^ 2 + (~b) ^ 2) / (~b) ^ 2) * Antiderivative(cos(~c + ~d * ~x) ^ (~m + 2) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n / sin(~(c') + ~(d') * ~x), ~x) => (-((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1)) / (~a * ~d * (~n + 1)) - (~b / (~a) ^ 2) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1), ~x)) + (1 / (~a) ^ 2) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 2) / sin(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n / cos(~(c') + ~(d') * ~x), ~x) => ((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1) / (~b * ~d * (~n + 1)) - (~a / (~b) ^ 2) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1), ~x)) + (1 / (~b) ^ 2) * Antiderivative((~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 2) / cos(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((((~a) ^ 2 + (~b) ^ 2) / (~a) ^ 2) * Antiderivative(sin(~c + ~d * ~x) ^ (~m + 2) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n, ~x) - ((2 * ~b) / (~a) ^ 2) * Antiderivative(sin(~c + ~d * ~x) ^ (~m + 1) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1), ~x)) + (1 / (~a) ^ 2) * Antiderivative(sin(~c + ~d * ~x) ^ ~m * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~n, -1) && LtQ(~m, -1)))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~m * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((((~a) ^ 2 + (~b) ^ 2) / (~b) ^ 2) * Antiderivative(cos(~c + ~d * ~x) ^ (~m + 2) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~n, ~x) - ((2 * ~a) / (~b) ^ 2) * Antiderivative(cos(~c + ~d * ~x) ^ (~m + 1) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 1), ~x)) + (1 / (~b) ^ 2) * Antiderivative(cos(~c + ~d * ~x) ^ ~m * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~n, -1) && LtQ(~m, -1)))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * sin(~(c') + ~(d') * ~x) ^ ~(n') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(p'), ~x) => Antiderivative(ExpandTrig(cos(~c + ~d * ~x) ^ ~m * sin(~c + ~d * ~x) ^ ~n * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * sin(~(c') + ~(d') * ~x) ^ ~(n') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (~a) ^ ~p * (~b) ^ ~p * Antiderivative(cos(~c + ~d * ~x) ^ ~m * sin(~c + ~d * ~x) ^ ~n * (~b * cos(~c + ~d * ~x) + ~a * sin(~c + ~d * ~x)) ^ -(~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && ILtQ(~p, 0))
    @apply_utils Antiderivative((cos(~(c') + ~(d') * ~x) ^ ~(m') * sin(~(c') + ~(d') * ~x) ^ ~(n')) / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => ((~b / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ ~m * sin(~c + ~d * ~x) ^ (~n - 1), ~x) + (~a / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ (~m - 1) * sin(~c + ~d * ~x) ^ ~n, ~x)) - ((~a * ~b) / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((cos(~c + ~d * ~x) ^ (~m - 1) * sin(~c + ~d * ~x) ^ (~n - 1)) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (IGtQ(~m, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((cos(~(c') + ~(d') * ~x) ^ ~(m') * sin(~(c') + ~(d') * ~x) ^ ~(n')) / (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => Antiderivative(ExpandTrig((cos(~c + ~d * ~x) ^ ~m * sin(~c + ~d * ~x) ^ ~n) / (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && IntegersQ(~m, ~n)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * ~x) ^ ~(m') * sin(~(c') + ~(d') * ~x) ^ ~(n') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~b / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ ~m * sin(~c + ~d * ~x) ^ (~n - 1) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~p + 1), ~x) + (~a / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ (~m - 1) * sin(~c + ~d * ~x) ^ ~n * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ (~p + 1), ~x)) - ((~a * ~b) / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(cos(~c + ~d * ~x) ^ (~m - 1) * sin(~c + ~d * ~x) ^ (~n - 1) * (~a * cos(~c + ~d * ~x) + ~b * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (IGtQ(~m, 0) && (IGtQ(~n, 0) && ILtQ(~p, 0))))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.5 trig^m (a cos+b sin)^n.jl


text:
## 05.jl
# n-gram

function ngram(input, n)
  words = []
  for i in 1:(length(input)-n+1)
    push!(words, input[i:i+n-1])
  end
  return words
end

println(ngram("I am an NLPer", 2))
println(ngram(utf32("私はNLPerです。"), 2))
println(ngram(split("I am an NLPer"," "), 2))
println(ngram(["私","は","NLPer","です","。"], 2))


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

# hygiene.jl
# Make sure that our macros have good hygiene

module M
using LinearAlgebra
using Test
import JuMP

model = JuMP.Model()
sense = JuMP.MathOptInterface.MIN_SENSE
JuMP.@variable(model, x >= 0)
r = 3:5
JuMP.@variable(model, y[i = r] <= i)
JuMP.@variable(model, z[i = 1:2, j = 1:2], Symmetric)
@test z isa LinearAlgebra.Symmetric

JuMP.@constraint(model, x + sum(j * y[j] for j in r) <= 1)
JuMP.@constraint(model, sum(y[j] for j = r if j == 4) <= 1)
JuMP.@constraint(model, -1 <= x + y[3] <= 1)
JuMP.@constraints(model, begin
    x + sum(j * y[j] for j in r) <= 1
    sum(y[j] for j = r if j == 4) <= 1
end)

JuMP.@SDconstraint(model, [x x; -x x] ⪰ 0)

JuMP.@objective(model, sense, y[4])
JuMP.@objective(model, Min, x + sum(j * y[j] for j in r))
JuMP.@objective(model, Max, sum(y[j] for j = r if j == 4))

JuMP.@NLconstraint(model, y[3] == 1)
JuMP.@NLconstraint(model, x + sum(j * y[j] for j in r) <= 1)
JuMP.@NLconstraint(model, sum(y[j] for j = r if j == 4) <= 1)
JuMP.@NLconstraints(model, begin
    x + sum(j * y[j] for j in r) <= 1
    sum(y[j] for j = r if j == 4) <= 1
end)

JuMP.@NLobjective(model, sense, y[4])
JuMP.@NLobjective(model, Min, x + sum(j * y[j] for j in r))
JuMP.@NLobjective(model, Max, sum(y[j] for j = r if j == 4))

# TODO: Add tests for the content of the model.

model = JuMP.Model()
i = 10
j = 10
JuMP.@expression(model, ex[j = 2:3], sum(i for i in 1:j))
@test ex[2] == 3
@test ex[3] == 6
@test i == 10
@test j == 10

end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#cd("Dropbox/PhD/JonAustin/Model")
using DelimitedFiles

# Include my functions.....
include("func_Gnp.jl")

#p = 0.1 # Probability of forming simplex
N = 50 # Number of nodes
m = 4 # simplex size
nvar_max = 1000
Nexp = 5

filename=string("Outputs/jl_4simplex","_N",N,"_m",m,"_Nexp",Nexp,".txt")

# Fill simplices with probability p
# -------------

# Keep track of simplices

for nvar=1:nvar_max

    p = 1/(nvar*nvar)
    final_results = Dict()

    println("doing run\t",nvar,"\tout of \t",nvar_max)
    for exp=1:Nexp
        Quadruples = Dict()
        Triples = Dict()
        Pairs = Dict()

        for a=0:N-m
            for b=a+1:N-(m-1)
                for c=b+1:N-(m-2)
                    for d=c+1:N-(m-3)
                        simplex = [a,b,c,d]
                        random_float = rand(Float64)
                        if (random_float<p)
                            Quadruples[simplex]=1
                            Triples = fill_triples(Triples,simplex)
                            Pairs = fill_pairs(Pairs,simplex)
                        end
                    end
                end
            end
        end

        # Evaluate prevalence of patterns
        # -------------
        results = Dict()
        for a=0:N-m
            for b=a+1:N-(m-1)
                for c=b+1:N-(m-2)
                    for d=c+1:N-(m-3)
                        simplex = [a,b,c,d]
                        results = find_pattern4(Quadruples,Triples,Pairs,results,simplex)
                    end
                end
            end
        end

        results = normalize_results(results)
        results["p"] = p
        final_results = archive_results(final_results,results,Nexp)
    end
    write_results(filename,final_results)
    #println(results)
end

max_stars_repo_name: jonassjuul/team-formation

max_issues_repo_name: jonassjuul/team-formation

max_forks_repo_name: jonassjuul/team-formation


text:
using MacroTools: @forward

Base.@kwdef struct I18nData
    data::Dict{String, Union{I18nData, String}} = Dict()
end
@forward I18nData.data Base.haskey, Base.keys, Base.values, Base.get, Base.get!, Base.getindex, Base.setindex!

Base.@kwdef mutable struct I18nConfig
    current_language::String = get_system_language()
    fallback::Vector{String} = String[]
end

Base.@kwdef mutable struct I18nContext{T <: Function}
    data::I18nData = I18nData()
    srcpath::String = ""
    override_config::Union{Nothing, I18nConfig} = nothing
    fallback::T
end
get_config(ctx::I18nContext) = if isnothing(ctx.override_config)
    GlobalI18nConfig
else
    ctx.override_config
end

const GlobalI18nConfig = I18nConfig()

const I18nContexts = Dict{Module, I18nContext}()

const OnLanguageChange = Dict{I18nConfig, Set{Function}}()


text:
using StatsBase, DataFrames, Flux
include("2048_fn.jl")

using Flux
model = Flux.Chain(
  Dense(14*16+4, 128, relu),
  Dense(128, 64, relu),
  Dense(64, 32, relu),
  Dense(32, 1, relu));

loss(x,y) = sum((model(x) .- y).^2)
opt = ADAM(params(model))

@time simulate_one_game_and_train!(model, opt, 0.9)

init_grid = init_game()
@time (init_grid, fnl, seq, cart, tf, reward) = simulate_game!(init_grid, model, 0.5)

# score every move using untrained model
grid_try = copy(init_grid)
initial_model = [Tracker.data(model(grid_move_to_feature(sim_seq(grid_try, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]

# now fit the model once
x = generate_x_based_on_episode(init_grid, seq, cart, tf);
y = discounted_rewards(reward, 0.99)[1:end-1];
xy = [(x[:,j], y[j]) for j in 1:size(x)[2]]
Flux.@epochs 100 Flux.train!(loss, xy, opt)

updated_model = [Tracker.data(model(grid_move_to_feature(sim_seq(grid_try, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]
using Plots
plot(cumsum(reward))
plot!(updated_model)
plot!(y)


text:
# Autogenerated wrapper script for ColPack_jll for powerpc64le-linux-gnu-cxx03
export libcolpack

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ColPack")
JLLWrappers.@declare_library_product(libcolpack, "libColPack.so.0")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcolpack,
        "lib/libColPack.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Ifsome

export ifsome
"""
    3|>ifsome(x->x+1) # 4
"""
ifsome(f1::Function) = _->ifsome(_,f1)

"nothing|>ifsome(x->x+1) # nothing"
ifsome(v::Void, f1::Function)=nothing

ifsome(sm, f::Function) = f(sm)

"3|>ifsome(x->x+1, ()->:no) # 4"
ifsome(f1::Function, f2::Function) = _->ifsome(_,f1,f2)

"nothing|>ifsome(x->x+1, ()->:no) # :no"
ifsome(v::Void, f1::Function, f2::Function) = f2()

ifsome(sm, f1::Function, f2::Function) = f1(sm)

"3|>ifsome(x->x+1, :no) # 4"
ifsome(f1::Function, v2::Any) = _->ifsome(_,f1,v2)

"nothing|>ifsome(x->x+1, :no) # :no"
ifsome(v::Void, f1::Function, v2)=v2

ifsome(sm, f1::Function, v2)=f1(sm)

end # module

max_stars_repo_name: closescreen/Ifsome

max_issues_repo_name: closescreen/Ifsome

max_forks_repo_name: closescreen/Ifsome


text:
#
# # OnlinePCA: Thomas Breloff (Cointegrated Technologies)
#
#
#
# # solving for a dimension-reduced Y = XV', where X (n x d) is the original data, and Y (n x k) is projected
# # V (k x d) is a matrix where the columns are the first k eigenvectors of X'X/n (the covariance of X)
# # e (k x 1) is a vector of eigenvectors of the covariance matrix
#
# # We compute e and V by incrementally updating e and U, where U is a (k x d) matrix with the properties:
# #		Uᵢ = iᵗʰ column of U
# #		Vᵢ = iᵗʰ column of V
# #			 = iᵗʰ eigenvector of X-covariance
# #			 = Uᵢ / ‖Uᵢ‖
# #		eᵢ = iᵗʰ eigenvalue of X-covariance
# #		   = ‖Uᵢ‖
#
# # note: V is the normalized version of U
#
# type OnlinePCA{W<:Weighting} <: OnlineStat
#
# 	d::Int  # number of input vars
#   k::Int  # number of principal components
#   weighting::W
#   n::Int
#
# 	U::MatF  # (k x d)
# 	V::MatF  # (k x d)
# 	e::VecF	 # (k x 1)
#   xmeans::Means{W}
#  end
#
# function OnlinePCA(d::Integer, k::Integer, wgt::Weighting = default(Weighting))
# 	OnlinePCA(d, k, wgt, 0, zeros(k,d), zeros(k,d), zeros(k), Means(d, wgt))
# end
#
#
# function OnlinePCA(x::AVecF, k::Integer, wgt::Weighting = default(Weighting))
# 	o = OnlinePCA(length(x), k, wgt)
# 	update!(o, x)
# 	o
# end
#
# function OnlinePCA(X::AMatF, k::Integer, wgt::Weighting = default(Weighting))
# 	o = OnlinePCA(ncols(X), k, wgt)
# 	update!(o, X)
# 	o
# end
#
#
#
# #-----------------------------------------------------------------------# state
#
# statenames(o::OnlinePCA) = [:U, :V, :e, :xmeans, :nobs]
# state(o::OnlinePCA) = Any[copy(o.U), copy(o.V), copy(o.e), copy(mean(o.xmeans)), nobs(o)]
#
#
# #---------------------------------------------------------------------# update!
#
# # TODO: roughly based on Weng et al (2003): "Candid covariance-free incremental principal component analysis"
# # used https://github.com/kevinhughes27/pyIPCA/blob/master/pyIPCA/ccipca.py as a reference
#
#
# # TODO: optimize, potentially by using view(X, ...) instead of X[...]
#
# function update!(o::OnlinePCA, x::AVecF)
#
# 	x = center!(o.xmeans, x)
# 	λ = weight(o)
#
# 	@inbounds for i in 1:min(o.k, o.n)
#
# 		if o.e[i] == 0. # this should be more robust than checking i == o.n
#
# 			# initialize ith principal component
# 			Ui = x
# 			ei = norm(Ui)
# 			Vi = Ui / ei
#
# 		else
#
# 			# update the ith eigvec/eigval
# 			Ui = row(o.U, i)
# 			smooth!(Ui, x * (dot(x, Ui) / o.e[i]), λ)
# 			ei = norm(Ui)
# 			Vi = Ui / ei
#
# 			# subtract projection on ith PC
# 			x -= dot(x, Vi) * Vi
#
# 		end
#
# 		# store these updates
# 		row!(o.U, i, Ui)
# 		row!(o.V, i, Vi)
# 		o.e[i] = ei
# 	end
#
# 	o.n += 1
# 	nothing
# end
#
#
# # function update!(o::OnlinePCA, X::AMatF)
# # 	for i in 1:nrows(X)
# # 		update!(o, row(X,i))
# # 	end
# # end
#
#
#
# function Base.empty!(o::OnlinePCA)
# 	o.U = zeros(o.k, o.d)
# 	o.V = zeros(o.k, o.d)
# 	o.e = zeros(o.k)
# 	o.n = 0
# 	o.xmeans = Means(o.d, o.weighting)
# 	nothing
# end
#
# function Base.merge!(o1::OnlinePCA, o2::OnlinePCA)
# 	error("Merging undefined for PCA")
# end
#
#
# # returns a vector z = Vx
# StatsBase.predict(o::OnlinePCA, x::AVecF) = o.V * center(o.xmeans, x)
#
# function StatsBase.predict(o::OnlinePCA, X::AMatF)
# 	n = size(X,1)
# 	Z = zeros(n, o.k)
# 	for i in 1:n
# 		row!(Z, i, predict(o, row(X,i)))
# 	end
# 	Z
# end


text:
"""

General Framework for the data augmented Gaussian Processes

"""
module AugmentedGaussianProcesses

    export AbstractGP, GP, VGP, SVGP, VStP, MCGP, MOSVGP
    export Likelihood,  RegressionLikelihood, ClassificationLikelihood, MultiClassLikelihood
    export GaussianLikelihood, StudentTLikelihood, LaplaceLikelihood, HeteroscedasticLikelihood
    export LogisticLikelihood, BayesianSVM
    export SoftMaxLikelihood, LogisticSoftMaxLikelihood
    export PoissonLikelihood, NegBinomialLikelihood
    export Inference, Analytic, AnalyticVI, AnalyticSVI, GibbsSampling, MCIntegrationVI, MCIntegrationSVI, QuadratureVI, QuadratureSVI
    export NumericalVI, NumericalSVI
    export PriorMean, ZeroMean, ConstantMean, EmpiricalMean, AffineMean

    #Useful functions and module
    include("kmeans/KMeansModule.jl")
    include("functions/PGSampler.jl")
    include("functions/GIGSampler.jl")
    #include("functions/PerturbativeCorrection.jl")
    # include("functions/GPAnalysisTools.jl")
    # include("functions/IO_model.jl")
    #Custom modules
    using .KMeansModule
    using .PGSampler
    using .GIGSampler
    # using .PerturbativeCorrection
    # using .GPAnalysisTools
    # using .IO_model
    #General modules
    using LinearAlgebra, Random
    using GradDescent
    using KernelFunctions
    using ForwardDiff
    using PDMats
    using AdvancedHMC, MCMCChains
    using DataFrames
    using StatsBase, StatsFuns, SpecialFunctions, Distributions, FastGaussQuadrature
    using ProgressMeter
    #Exported modules
    export KMeansModule
    export KMeansInducingPoints
    #Useful functions
    export train!, sample
    export predict_f, predict_y, proba_y
    export fstar, ELBO
    export covariance, diag_covariance, prior_mean

    # Main classes
    abstract type Inference{T<:Real} end
    abstract type VariationalInference{T} <: Inference{T} end
    abstract type SamplingInference{T} <: Inference{T} end
    abstract type Likelihood{T<:Real} end
    abstract type Abstract_GP{T<:Real} end

    const LatentArray = Vector #For future optimization : How collection of latent GP parameters and local variables are stored
    include("prior/priormean.jl")
    include("models/inducing_points.jl")

    include("models/AbstractGP.jl")
    include("models/GP_base.jl")
    include("models/GP.jl")
    include("models/VGP.jl")
    include("models/MCGP.jl")
    include("models/SVGP.jl")
    include("models/VStP.jl")
    include("models/MOSVGP.jl")

    include("inference/inference.jl")
    include("likelihood/likelihood.jl")

    include("functions/utils.jl")
    include("functions/init.jl")
    include("functions/KLdivergences.jl")
    include("functions/plotting.jl")

    #Training Functions
    include("training.jl")
    include("autotuning.jl")
    include("predictions.jl")
end #End Module