text:
function setupJointInversion(m,filenamePrefix::String,plotting::Bool,
		workersFWI::Array{Int64,1}=workers(),maxBatchSize::Int64 = 48,
		Ainv::AbstractSolver = getMUMPSsolver([],0,0,2), misfun::Function=SSDFun,useFilesForFields::Bool = false,onlyFWI::Bool=true,beta::Float64=1.0)

Q,P,pMisFWIRFs,SourcesSubIndFWI,contDivFWI,Iact,sback,mref,boundsHigh,boundsLow = setupFWI(m,filenamePrefix,plotting,workersFWI,maxBatchSize,Ainv, misfun,useFilesForFields)

if onlyFWI==false
	file = matopen(string(filenamePrefix,"_PARAM.mat"));
	n_cells = read(file,"n");
	OmegaDomain = read(file,"domain");
	Minv = getRegularMesh(OmegaDomain,n_cells);
	HO = read(file,"HO");
	close(file);
	### Read receivers and sources files
	RCVfile = string(filenamePrefix,"_rcvMap.dat");
	SRCfile = string(filenamePrefix,"_srcMap.dat");
	srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
	rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);

	println("Reading Eik data:");

	(DobsEik,WdEik) =  readDataFileToDataMat(string(filenamePrefix,"_travelTime.dat"),srcNodeMap,rcvNodeMap);

	EikMPIWorkers = nworkers(); # this just set the maximal MPI workers. To activate parallelism, run addprocs()

	(pForEik,contDivEik,SourcesSubIndEik) = getEikonalInvParam(Minv,Q,P,HO,EikMPIWorkers,useFilesForFields);
	Wd 		= Array{Array{Float64}}(length(pForEik));
	dobs 	= Array{Array{Float64}}(length(pForEik));
	for i=1:length(pForEik)
		I_i = SourcesSubIndEik[i];
		Wd[i]   = beta*WdEik[:,I_i];
		dobs[i] = DobsEik[:,I_i];
	end
	WdEik = 0;
	# DobsEik = 0;

	pMisEikRFs = getMisfitParam(pForEik, Wd, dobs, misfun, Iact,sback);
end

# ########################################################################################################
# ##### Set up remote workers ############################################################################
# ########################################################################################################

if onlyFWI==false
	SourcesSubInd    = [SourcesSubIndEik;SourcesSubIndFWI]; SourcesSubIndEik = 0;SourcesSubIndFWI = 0;
	pMis 			 = [pMisEikRFs;pMisFWIRFs];pMisEikRFs = 0;pMisFWIRFs = 0;
	contDiv			 = [contDivEik[1];contDivFWI + contDivEik[2] - 1];contDivEik = 0;contDivFWI = 0;
else
	SourcesSubInd    = SourcesSubIndFWI;
	pMis 			 = pMisFWIRFs;
	contDiv			 = contDivFWI;
end

return Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow;

end

# function setupJointInversion(m,filenamePrefix::String,plotting::Bool,
		# workersFWI::Array{Int64,1}=workers(),maxBatchSize::Int64 = 48,
		# Ainv::AbstractSolver = getMUMPSsolver([],0,0,2),onlyFWI::Bool=true, misfun::Function=SSDFun,beta::Float64=1.0,useFilesForFields::Bool = false)

# file = matopen(string(filenamePrefix,"_PARAM.mat"));
# n_cells = read(file,"n");
# OmegaDomain = read(file,"domain");
# Minv = getRegularMesh(OmegaDomain,n_cells);
# gamma = read(file,"gamma");
# omega = read(file,"omega");
# waveCoef = read(file,"waveCoef");
# if length(omega)==1
	# # matlab saves a 1 variable array as scalar.
	# omega = [omega];
	# waveCoef = [waveCoef];
# end
# println("Note: Assuming wavelet coeficient is taken into account in the observed data. RHS in FWI is just a delta func.");
# waveCoef = ones(ComplexF64,size(waveCoef))

# if onlyFWI==false
	# HO = read(file,"HO");
# end

# boundsLow = read(file,"boundsLow");
# boundsHigh = read(file,"boundsHigh");
# mref =  read(file,"mref");
# close(file);

# ### Read receivers and sources files
# RCVfile = string(filenamePrefix,"_rcvMap.dat");
# SRCfile = string(filenamePrefix,"_srcMap.dat");

# srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
# rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);

# Q = generateSrcRcvProjOperators(Minv.n.+1,srcNodeMap);
# Q = Q.*1/(norm(Minv.h)^2);
# println("We have ",size(Q,2)," sources");
# P = generateSrcRcvProjOperators(Minv.n.+1,rcvNodeMap);

# ########################################################################################################
# ##### Set up remote workers ############################################################################
# ########################################################################################################

# N = prod(Minv.n);

# Iact = SparseMatrixCSC(1.0I,N,N);
# mback   = zeros(Float64,N);
# ## Setting the sea constant:
# mask = zeros(N);
# mask[abs.(m[:] .- minimum(m)) .< 1e-2] .= 1.0;
# mask[gamma[:] .>= 0.95*maximum(gamma)] .= 1.0;
# # setup active cells
# mback = vec(mref[:].*mask);
# # mback = vec(m[:].*mask);
# sback = velocityToSlowSquared(mback)[1];
# sback[mask .== 0.0] .= 0.0;
# Iact = Iact[:,mask .== 0.0];

# boundsLow = Iact'*boundsLow;
# boundsHigh = Iact'*boundsHigh;
# mref = Iact'*mref[:];

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

# if onlyFWI==false
	# println("Reading Eik data:");
	# (DobsEik,WdEik) =  readDataFileToDataMat(string(filenamePrefix,"_travelTime.dat"),srcNodeMap,rcvNodeMap);

	# EikMPIWorkers = nworkers(); # this just set the maximal MPI workers. To activate parallelism, run addprocs()

	# (pForEik,contDivEik,SourcesSubIndEik) = getEikonalInvParam(Minv,Q,P,HO,EikMPIWorkers,useFilesForFields);
	# Wd 		= Array{Array{Float64}}(length(pForEik));
	# dobs 	= Array{Array{Float64}}(length(pForEik));
	# for i=1:length(pForEik)
		# I_i = SourcesSubIndEik[i];
		# Wd[i]   = beta*WdEik[:,I_i];
		# dobs[i] = DobsEik[:,I_i];
	# end
	# WdEik = 0;
	# # DobsEik = 0;

	# pMisEikRFs = getMisfitParam(pForEik, Wd, dobs, misfun, Iact,sback);
# end
# println("Reading FWI data:");
# batch = min(size(Q,2),maxBatchSize);
# (pForFWI,contDivFWI,SourcesSubIndFWI) = getFWIparam(omega,waveCoef,vec(gamma),Q,P,Minv,Ainv,workersFWI,batch,useFilesForFields);

# # write data to remote workers
# Wd   = Array{Array{ComplexF64,2}}(undef,length(pForFWI))
# dobs = Array{Array{ComplexF64,2}}(undef,length(pForFWI))
# for k = 1:length(omega)
	# omRound = string(round((omega[k]/(2*pi))*100.0)/100.0);
	# (DobsFWIwk,WdFWIwk) =  readDataFileToDataMat(string(filenamePrefix,"_freq",omRound,".dat"),srcNodeMap,rcvNodeMap);

	# for i = contDivFWI[k]:contDivFWI[k+1]-1
		# I_i = SourcesSubIndFWI[i]; # subset of sources for ith worker.
		# Wd[i] 	= WdFWIwk[:,I_i];
		# # Wd[i]   = exp(-1im*omega[k]*DobsEik[:,I_i])
		# dobs[i] = DobsFWIwk[:,I_i];

	# end
	# DobsFWIwk = 0;
	# WdFWIwk = 0;
# end

# pMisFWIRFs = getMisfitParam(pForFWI, Wd, dobs, misfun, Iact,sback);

# ########################################################################################################
# ##### Set up remote workers ############################################################################
# ########################################################################################################

# if onlyFWI==false
	# SourcesSubInd    = [SourcesSubIndEik;SourcesSubIndFWI]; SourcesSubIndEik = 0;SourcesSubIndFWI = 0;
	# pMis 			 = [pMisEikRFs;pMisFWIRFs];pMisEikRFs = 0;pMisFWIRFs = 0;
	# contDiv			 = [contDivEik[1];contDivFWI + contDivEik[2] - 1];contDivEik = 0;contDivFWI = 0;
# else
	# SourcesSubInd    = SourcesSubIndFWI;
	# pMis 			 = pMisFWIRFs;
	# contDiv			 = contDivFWI;
# end

# return Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow;
# end


text:
import ClimateMachine.BalanceLaws:
    # declaration
    vars_state,
    # initialization
    nodal_init_state_auxiliary!,
    init_state_prognostic!,
    init_state_auxiliary!,
    # rhs computation
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
    # boundary conditions
    boundary_conditions,
    boundary_state!
import ClimateMachine.NumericalFluxes: numerical_boundary_flux_first_order!

struct DryReferenceState{TP}
    temperature_profile::TP
end

"""
    Declaration of state variables
    vars_state returns a NamedTuple of data types.
"""
function vars_state(m::DryAtmosModel, st::Auxiliary, FT)
    @vars begin
        x::FT
        y::FT
        z::FT
        Φ::FT
        ∇Φ::SVector{3, FT} # TODO: only needed for the linear model
        ref_state::vars_state(m, m.physics.ref_state, st, FT)
        T_sfc::FT  # Ocean surface temperature
    end
end

vars_state(::DryAtmosModel, ::DryReferenceState, ::Auxiliary, FT) =
    @vars(T::FT, p::FT, ρ::FT, ρu::SVector{3, FT}, ρe::FT, ρq::FT)
vars_state(::DryAtmosModel, ::NoReferenceState, ::Auxiliary, FT) = @vars()

function vars_state(::DryAtmosModel, ::Prognostic, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        ρq::FT
        F_ρe_accum::FT # accumulated energy flux
    end
end

"""
    Initialization of state variables
    init_state_xyz! sets up the initial fields within our state variables
    (e.g., prognostic, auxiliary, etc.), however it seems to not initialized
    the gradient flux variables by default.
"""
function init_state_prognostic!(model::DryAtmosModel, state::Vars, aux::Vars, localgeo, t)
    x = aux.x
    y = aux.y
    z = aux.z

    parameters = model.physics.parameters
    ic = model.initial_conditions

    # TODO!: Set to 0 by default or assign IC
    if !isnothing(ic)
        state.ρ = ic.ρ(parameters, x, y, z)
        state.ρu = ic.ρu(parameters, x, y, z)
        state.ρe = ic.ρe(parameters, x, y, z)
        state.ρq = ic.ρq(parameters, x, y, z)
    end

    state.F_ρe_accum = 0

    return nothing
end

function nodal_init_state_auxiliary!(model::DryAtmosModel, state_auxiliary, tmp, geom)
    init_state_auxiliary!(model, model.physics.orientation, state_auxiliary, geom)
    init_state_auxiliary!(model, model.physics.ref_state, state_auxiliary, geom)

    FT = eltype(state_auxiliary)
    state_auxiliary.T_sfc = FT(0) # otherwise can get NaNs if only the boundary of the MPIStateArray is allocated
end

function init_state_auxiliary!(model::DryAtmosModel, ::SphericalOrientation, state_auxiliary, geom)
    g = model.physics.parameters.g

    r = norm(geom.coord)
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]
    state_auxiliary.Φ = g * r
    state_auxiliary.∇Φ = g * geom.coord / r
end

function init_state_auxiliary!(model::DryAtmosModel, ::FlatOrientation, state_auxiliary, geom)
    g = model.physics.parameters.g

    FT = eltype(state_auxiliary)

    r = geom.coord[3]
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]
    state_auxiliary.Φ = g * r
    state_auxiliary.∇Φ = SVector{3, FT}(0, 0, g)
end

function init_state_auxiliary!(::DryAtmosModel, ::NoReferenceState, state_auxiliary, geom) end

function init_state_auxiliary!(model::DryAtmosModel, ref_state::DryReferenceState, state_auxiliary, geom)
    orientation = model.physics.orientation
    R_d = model.physics.parameters.R_d
    γ = model.physics.parameters.γ
    Φ = state_auxiliary.Φ

    FT = eltype(state_auxiliary)

    # Calculation of a dry reference state
    z = altitude(model, orientation, geom)
    T, p = ref_state.temperature_profile(model.physics.parameters, z)
    ρ = p / R_d / T
    ρu = SVector{3, FT}(0, 0, 0)
    ρe = p / (γ - 1) + dot(ρu, ρu) / 2ρ + ρ * Φ
    ρq = FT(0)

    state_auxiliary.ref_state.T = T
    state_auxiliary.ref_state.p = p
    state_auxiliary.ref_state.ρ = ρ
    state_auxiliary.ref_state.ρu = ρu
    state_auxiliary.ref_state.ρe = ρe
    state_auxiliary.ref_state.ρq = ρq
end

"""
    LHS computations
"""
@inline function flux_first_order!(model::DryAtmosModel, flux::Grad, state::Vars, aux::Vars, t::Real, direction)
    lhs = model.physics.lhs
    physics = model.physics

    ntuple(Val(length(lhs))) do s
        Base.@_inline_meta
        calc_component!(flux, lhs[s], state, aux, physics)
    end
end

"""
    RHS computations
"""
function source!(model::DryAtmosModel, source, state_prognostic, state_auxiliary, _...)
    sources = model.physics.sources
    physics = model.physics

    ntuple(Val(length(sources))) do s
        Base.@_inline_meta
        calc_component!(source, sources[s], state_prognostic, state_auxiliary, physics)
    end

end

"""
    Boundary conditions with defaults
"""
boundary_conditions(model::DryAtmosModel) = model.boundary_conditions

function boundary_state!(::CentralNumericalFluxSecondOrder, _...)
    nothing
end

"""
    Utils
"""
function altitude(model::DryAtmosModel, ::SphericalOrientation, geom)
    return norm(geom.coord) - model.physics.parameters.a
end

function altitude(::DryAtmosModel, ::FlatOrientation, geom)
    @inbounds geom.coord[3]
end

ext: jl

max_stars_repo_name: CliMA/ClimaCoupler


text:
function diagnose_u(ψ, d::Domain)
    u = new_field(d)
    for I = CartesianIndices(d)
        u[I] = -∂y(I,ψ,d)
    end
    return u
end

function diagnose_v(ψ, d::Domain)
    v = new_field(d)
    for I = CartesianIndices(d)
        v[I] = ∂x(I,ψ,d)
    end
    return v
end

function diagnose_umag(ψ, d::Domain)
    umag = new_field(d)
    for I = CartesianIndices(d)
        umag[I] = sqrt(∂y(I,ψ,d)^2 + ∂x(I,ψ,d)^2)
    end
    return umag
end

function diagnose_ζ(ψ, d::Domain)
    ζ = new_field(d)
    for I = CartesianIndices(d)
        ζ[I] = ∂x²(I,ψ,d) + ∂y²(I,ψ,d)
    end
    return ζ
end

function diagnose_θ(ϕ, d::Domain)
    θ = new_field(d)
    for I = CartesianIndices(d)
        θ[I] = -∂z(I,ϕ,d)
    end
    return θ
end

function diagnose_penetration_depth(ψ, q, d::Domain)
    qmax, Imax = findmax(q)
    kmax = Imax[3]
    umag = diagnose_umag(ψ, d)
    ubot = umag[:,:,1]
    umax = umag[:,:,kmax]
    return maximum(ubot)/maximum(umax)
end


text:
using Flux
println("Integro Differential equation tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
println("Starting Soon!")
using GalacticOptim
using Optim
using Quadrature,Cubature, Cuba
using QuasiMonteCarlo
using SciMLBase
import ModelingToolkit: Interval, infimum, supremum
using DomainSets

using Random
Random.seed!(100)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

#Integration Tests
println("Integral Tests")
@parameters t
@variables i(..)
Di = Differential(t)
Ii = Integral(t in DomainSets.ClosedInterval(0, t))
eq = Di(i(t)) + 2*i(t) + 5*Ii(i(t)) ~ 1
bcs = [i(0.) ~ 0.0]
domains = [t ∈ Interval(0.0,2.0)]
chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[t],[i(t)])
sym_prob = NeuralPDE.symbolic_discretize(pde_system, discretization)
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
ts = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]
phi = discretization.phi

analytic_sol_func(t) = 1/2*(exp(-t))*(sin(2*t))
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(phi([t],res.minimizer)) for t in ts]
@test Flux.mse(u_real, u_predict) < 0.001
# plot(ts,u_real)
# plot!(ts,u_predict)

## Simple Integral Test
println("Simple Integral Test")

@parameters x
@variables u(..)
Ix = Integral(x in DomainSets.ClosedInterval(0, x))
# eq = Ix(u(x)) ~ (x^3)/3
eq = Ix(u(x)*cos(x))~ (x^3)/3

bcs = [u(0.) ~ 0.0]
domains = [x ∈ Interval(0.0,1.00)]
# chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))
chain = FastChain(FastDense(1,15,Flux.σ),FastDense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = initθ,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[x],[u(x)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=200)
xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]
phi = discretization.phi
u_predict  = [first(phi([x],res.minimizer)) for x in xs]
u_real  = [x^2/cos(x) for x in xs]
@test Flux.mse(u_real, u_predict) < 0.001

# plot(xs,u_real)
# plot!(xs,u_predict)

#simple multidimensitonal integral test
println("simple multidimensitonal integral test")

@parameters x,y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Ix = Integral((x,y) in DomainSets.UnitSquare())
eq = Ix(u(x,y)) ~ 1/3
bcs = [u(0., 0.) ~ 1, Dx(u(x,y)) ~ -2*x , Dy(u(x ,y)) ~ -2*y ]
domains = [x ∈ Interval(0.0,1.00), y ∈ Interval(0.0,1.00)]
chain = Chain(Dense(2,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
xs = 0.00:0.01:1.00
ys = 0.00:0.01:1.00
phi = discretization.phi

u_real = collect(1 - x^2 - y^2 for y in ys, x in xs);
u_predict = collect(Array(phi([x,y], res.minimizer))[1] for y in ys, x in xs);
@test Flux.mse(u_real, u_predict) < 0.001

# error_ = u_predict .- u_real
# p1 = plot(xs,ys,u_real,linetype=:contourf,label = "analytic")
# p2 = plot(xs,ys,u_predict,linetype=:contourf,label = "predict")
# p3 = plot(xs,ys,error_,linetype=:contourf,label = "error")
# plot(p1,p2,p3)

@parameters x,y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Ix = Integral((x,y) in DomainSets.ProductDomain(UnitInterval(),ClosedInterval(0 ,x)))
eq = Ix(u(x,y)) ~ 5/12
bcs = [u(0., 0.) ~ 0, Dy(u(x,y)) ~ 2*y , u(x, 0) ~ x ]
domains = [x ∈ Interval(0.0,1.00), y ∈ Interval(0.0,1.00)]
chain = Chain(Dense(2,15,Flux.σ),Dense(15,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy_;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             )
@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)
xs = 0.00:0.01:1.00
ys = 0.00:0.01:1.00
phi = discretization.phi

u_real = collect( x + y^2 for y in ys, x in xs);
u_predict = collect(Array(phi([x,y], res.minimizer))[1] for y in ys, x in xs);
@test Flux.mse(u_real, u_predict) < 0.01

# error_ = u_predict .- u_real
# p1 = plot(xs,ys,u_real,linetype=:contourf,label = "analytic")
# p2 = plot(xs,ys,u_predict,linetype=:contourf,label = "predict")
# p3 = plot(xs,ys,error_,linetype=:contourf,label = "error")
# plot(p1,p2,p3)

## Two variables Integral Test
println("Two variables Integral Test")

@parameters x
@variables u(..) w(..)
Dx = Differential(x)
Ix = Integral(x in DomainSets.ClosedInterval(1, x))
# Iinf = Integral(x in DomainSets.ClosedInterval(1, Inf))

eqs =  [Ix(u(x)*w(x)) ~ log(abs(x)),
        Dx(w(x)) ~ -2/(x^3),
        u(x) ~ x ]

bcs = [u(1.) ~ 1.0, w(1.) ~ 1.0]
domains = [x ∈ Interval(1.0,2.0)]

chains = [FastChain(FastDense(1,15,Flux.σ),FastDense(15,1)) for _ in 1:2]
initθ = map(chain -> Float64.(DiffEqFlux.initial_params(chain)),chains)
strategy_ = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chains,
                                             strategy_;
                                             init_params = initθ
                                             )
@named pde_system = PDESystem(eqs,bcs,domains,[x],[u(x), w(x)])
prob = NeuralPDE.discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=200)
xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]
phi = discretization.phi
initθ = discretization.init_params
acum =  [0;accumulate(+, length.(initθ))]
sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
minimizers = [res.minimizer[s] for s in sep]

u_predict  = [(phi[1]([x],minimizers[1]))[1] for x in xs]
w_predict  = [(phi[2]([x],minimizers[2]))[1] for x in xs]
u_real  = [x for x in xs]
w_real  = [1/x^2 for x in xs]
@test Flux.mse(u_real, u_predict) < 0.001
@test Flux.mse(w_real, w_predict) < 0.001

# plot(xs,u_real)
# plot!(xs,u_predict)
# plot(xs,w_real)
# plot!(xs,w_predict)


text:
module STBUnits
    #
    # Not Used presently - initial attempt to use Unitful to model all the values in hhlds and parameters.
    #
    #
    using Unitful

    @unit money "GBP" Money 1u"£" false
    @unit pA "PA" PA 365u"£" false;
    @unit rate "r" Rate 1u"r" false

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
@testset "Plotting" begin
  Random.seed!(123)
  z₁ = randn(10000)
  z₂ = z₁ + randn(10000)
  d = georef(DataFrame(z₁=z₁,z₂=z₂), RegularGrid(100,100))

  @testset "DistPlot1D" begin
    if visualtests
      @plottest distplot1d(d,:z₁) joinpath(datadir,"distplot1D.png") !isCI
    end
  end

  @testset "DistPlot2D" begin
    if visualtests
      @plottest distplot2d(d,:z₁,:z₂) joinpath(datadir,"distplot2D.png") !isCI
    end
  end

  @testset "CornerPlot" begin
    if visualtests
      @plottest cornerplot(d) joinpath(datadir,"cornerplot.png") !isCI
    end
  end

  @testset "HScatter" begin
    if visualtests
      @plottest begin
        sdata = readgeotable(joinpath(datadir,"samples2D.tsv"), coordnames=(:x,:y))
        p0 = hscatter(sdata, :value, lag=0)
        p1 = hscatter(sdata, :value, lag=1)
        p2 = hscatter(sdata, :value, lag=2)
        p3 = hscatter(sdata, :value, lag=3)
        plot(p0, p1, p2, p3, layout=(2,2), size=(600,600))
      end joinpath(datadir,"hscatter.png") !isCI
    end
  end
end


text:
#
# Functions used to reduce the lists in parallel calculations
#

"""
```
keep_best_list!(list1,list2)
```

Update `list1` considering the list of minimum distances given in `list2`.

Internal function or structure - interface may change.

"""
function keep_best_list!(list1, list2)
    for i in eachindex(list1)
        if list2[i].d < list1[i].d
            list1[i] = list2[i]
        end
    end
end

"""

```
reduce_list!(list, list_threaded)
```

Update the final `list` of minimum-distances given the threaded list `list_threaded`.

Internal function or structure - interface may change.

"""
function reduce_list!(list, list_threaded)
    list .= list_threaded[1]
    for it in 2:length(list_threaded)
        keep_best_list!(list, list_threaded[it])
    end
    return list
end

"""

```
reduce_list_pair!(lists, lists_threaded)
```

Update the final tuple `lists`, of minimum-distances of the two molecule sets given the threaded lists `lists_threaded`.

Internal function or structure - interface may change.

"""
function reduce_list_pair!(lists, lists_threaded)
    lists[1] .= lists_threaded[1][1]
    lists[2] .= lists_threaded[1][2]
    for it in 2:length(lists_threaded)
        keep_best_list!(lists[1], lists_threaded[it][1])
        keep_best_list!(lists[2], lists_threaded[it][2])
    end
    return lists
end


text:
#@testset "BigMat.jl" begin
    using LinearAlgebra
    # tests for \ (solving of linear system Ax = b and A*inv(K) = I)

    A = rand(64,64);
    A = A*A'/2;
    b = rand(64,1);

    mps,err = MPT_SVD(b,[4 4 4],0.0);
    mpo,err = MPT_SVD(A,[4 4 4; 4 4 4],0.0)
    mpo\mps
    @test norm(mps2vec(mpo\mps)-A\b)/norm(A\b) < 1e-5

    # solve directly for inverse with A*inv(A) = I
    id        = eye([4 4 4; 4 4 4]);
    mpo1      = MPT([mpo[1:3]..., id[1:3]...]);
    id        = Matrix(I,64,64);
    id        = 1.0*id[:];
    idvec,err = MPT_SVD(id,[4 4 4 4 4 4],0.0);
    invA      = mpo1\idvec;
    reshape(mps2vec(invA),(64,64));
    @test norm(lyap-inv(A))/norm(inv(A)) < 1e-5

    # solve Lyaponov equation ((I ⊗ A)+(A^T ⊗ I))vec(X) = 2vec(I)
    id = eye([4 4 4; 4 4 4]);
    mpo1 = MPT([mpo[1:3]..., id[1:3]...]);
    mpoT = transpose(mpo);
    mpo2 = MPT([id[1:3]..., mpoT[1:3]...]);

    sum = mpo1 + mpo2;
    id = Matrix(I,64,64);
    id = 2.0*id[:];
    idvec,err = MPT_SVD(id,[4 4 4 4 4 4],0.0);

    lyap = sum\idvec;
    lyap = reshape(mps2vec(lyap),(64,64));
    @test norm(lyap-inv(A))/norm(inv(A)) < 1e-5

    # test of column-wise and row-wise Kronecker product (KathriRao)
    A = rand(3,4);
    B = rand(3,4);
    C = rand(3,4);
    @test size(KathriRao(A,B,1),2) == 16
    @test size(KathriRao(A,B,2),1) == 9

    @test norm(A.*B.*C) < norm(KathriRao(A,KathriRao(B,C,2),2))

    # test matrix multiplication
    @test norm(mps2vec(mpo*mpo*mps) - A*A*b)/norm(A*A*b) < 1e-10

    A = randn(8,8)
    B = randn(8,8)
    Attm,err = MPT_SVD(A,[2 2 2;2 2 2],0.0)
    AB = *(Attm,B,3)
    norm(AB-A*B)/norm(A*B)
    ABA = AB*transpose(Attm)
    norm(ABA-A*B*A')/norm(A*B*A')

    # test pseudo inverse
    @run approxpseudoinverse(Attm,0.0,0.0)
    P = approxpseudoinverse(Attm,0.0,0.0)
    norm(inv(A) - mpo2mat(P))/norm(inv(A))

    # test vectorbymatrix
    A = randn(8,2)
    Attm,err = MPT_SVD(A,[2 2 2;1 1 2],0.0)
    b = randn(8)
    btt,err = MPT_SVD(b,[2 2 2],0.0)
    norm(vectorbymatrix(Attm,btt) - b'*A)/norm(b'*A)

    # test matrixbyvector
    matrixbyvector(Attm,btt)
    norm(matrixbyvector(Attm,btt) - A'*b)/norm(A'*b)

#end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using Random

#TODO for some reason Sarsa did not work with this policy. Used a q function and mual eps greed choice
mutable struct EpsilonGreedy{TP, TB} <: AbstractPolicy where {TP<:Real,TB<:AbstractFuncApprox}
    f::LinearFunction{TP, TB}
    action::Int
    vals::Array{TP,1}
    ϵ::TP
    function EpsilonGreedy(state_dim::Int, num_actions::Int, ϵ::T) where {T<:Real}
        f = LinearFunction(state_dim, num_actions)
        new{eltype(f.θ),typeof(f.ϕ)}(f, -1, zeros(eltype(f.θ), num_actions))
    end
    function EpsilonGreedy(::Type{T}, state_dim::Int, num_actions::Int, ϵ::T) where {T<:Real}
        f = LinearFunction(T, state_dim, num_actions)
        new{T,typeof(f.ϕ)}(f, -1, zeros(T, num_actions), ϵ)
    end
    function EpsilonGreedy(f::LinearFunction{T}, ϵ::T) where {T<:Real}
        new{T,typeof(f.ϕ)}(clone(f), -1, zeros(T, get_num_outputs(f)), ϵ)
    end
end

function get_num_params(π::EpsilonGreedy)::Int
    return get_num_params(π.f)
end

function get_action!(π::EpsilonGreedy{T}, x, rng::AbstractRNG)::T where {T<:Real}
    call!(π.vals, π.f, x)
    num_actions = length(π.vals)
    logp = π.ϵ * (1. / T(num_actions))
    if rand(rng, 1)[1] < π.ϵ
        π.action = argmax(π.vals)
        logp += (1. - π.ϵ)
    else
        π.action = rand(rng, 1:length(π.vals))[1]
    end
    return logp
end

function gradient_logp!(grad::Array{T}, π::EpsilonGreedy{T}, x, action::Int)::T where {T<:Real}
    grad = reshape(grad, size(get_params(π.f)))
    gradient!(grad, π.f, x, action)
    logp = 0.
    return logp
end

function get_action_gradient_logp!(grad::Array{T}, π::EpsilonGreedy{T}, x, rng::AbstractRNG)::T where {T<:Real}
    call_gradient!(π.vals, grad, π.f, x)
    num_actions = length(π.vals)
    logp = π.ϵ * (1. / T(num_actions))
    if rand(rng, 1)[1] < π.ϵ
        π.action = argmax(π.vals)
        logp += (1. - π.ϵ)
    else
        π.action = rand(rng, 1:length(π.vals))[1]
    end

    fill!(view(grad, :, 1:num_actions .!= π.action), 0.)
    return logp
end

function get_action_gradient_logp(grad::Array{T}, π::EpsilonGreedy{T}, x, rng::AbstractRNG) where {T<:Real}
    grad = zeros(T, size(get_params(π)))
    logp = get_action_gradient_logp!(grad, π, x)
    action = π.action
    return action, grad, logp
end

function set_params!(π::EpsilonGreedy{T}, θ::Array{T}) where {T<:Real}
    set_params!(π.f, θ)
end

function get_params(π::EpsilonGreedy)
    get_params(π.f)
end

function copy_params!(params::Array{T}, π::EpsilonGreedy{T}) where {T<:Real}
    copy_params!(params, π.f)
end

function copy_params(π::EpsilonGreedy{T})::Array{T} where {T<:Real}
    return copy_params(π.f)
end

function add_to_params!(π::EpsilonGreedy{T}, Δθ::Array{T}) where {T<:Real}
    add_to_params!(π.f, Δθ)
end

function clone(π::EpsilonGreedy)::EpsilonGreedy
    π₂ = EpsilonGreedy(π.f, π.ϵ)
    π₂.action = deepcopy(π.action)
    π₂.vals = deepcopy(π.vals)
    return π₂
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# LibGit2.jl "clone" example - shows how to use the clone API
#
# Written by the LibGit2.jl contributors, based on the libgit2 examples.
#
# To the extent possible under law, the author(s) have dedicated all copyright
# and related and neighboring rights to this software to the public domain
# worldwide. This software is distributed without any warranty.
#
# You should have received a copy of the CC0 Public Domain Dedication along
# with this software. If not, see
# <http://creativecommons.org/publicdomain/zero/1.0/>.

using LibGit2
using ArgParse

# This example demonstrates the use of the LibGit2.jl clone APIs
# emulating the command `git clone` and printing the progress of
# the clone process and of the checkout

# from the git-clone man page:
# Clones a repository into a newly created directory, creates remote-tracking
# branches for each branch in the cloned repository (visible using git branch -r),
# and creates and checks out an initial branch that is forked from the cloned
# repository's currently active branch.

function parse_commandline()
    settings = ArgParseSettings(autofix_names=true)
    @add_arg_table settings begin
        "repo-url"
        help = "repo url"
        required = true

        "directory"
        help = ""
    end
    return parse_args(settings)
end

# Assuming n is in bytes
function pretty_size(n)
    c = log2(n) / 10
    val = n / 1024^floor(c)
    c < 1      && return @sprintf("%i B",val)
    1 <= c < 2 && return @sprintf("%.2f KiB",val)
    2 <= c < 3 && return @sprintf("%.2f MiB",val)
    3 <= c     && return @sprintf("%.2f GiB",val)
end

# fetch and checkout head with progress indication
function doclone(repo_url::AbstractString, dest_dir::AbstractString)
    if isempty(dest_dir)
        git_dir = joinpath(pwd(), basename(splitext(repo_url)[1]))
    else
        git_dir = joinpath(pwd(), basename(dest_dir))
    end

    info("cloning into $git_dir")
    mkdir(git_dir)

    total_objects = indexed_objects = received_objects = received_bytes = 0
    repo = LibGit2.repo_clone(repo_url, git_dir, Dict{Any,Any}(
        :callbacks => Dict{Any,Any}(
          :transfer_progress => (args...) -> begin
                total_objects, indexed_objects, received_objects, received_bytes = args
                perc = int(received_objects / total_objects * 100.0)
                print("\rReceiving objects: $perc% ($received_objects/$total_objects) ",
                      "$(pretty_size(received_bytes))")
            end
          )
        ))

    println(", done.")

    path = completed_steps = total_steps = payload = 0
    LibGit2.checkout_head!(repo,
        Dict{Any,Any}(
         :strategy => :safe_create,
         :progress => (args...) -> begin
            path, completed_steps, total_steps = args
            print("\rUnpacking files: $completed_steps/$total_steps")
         end
        ))

    println(", done.")
    return repo
end

function main()
    opts = parse_commandline()
    if haskey(opts, "directory") && opts["directory"] != nothing
        doclone(opts["repo_url"], opts["directory"])
    else
        doclone(opts["repo_url"], "")
    end
end

main()


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

__precompile__(true)

module Unicode

using Base.Unicode: normalize, graphemes, isassigned, iscased

export graphemes

# BEGIN 0.7 deprecations

@deprecate is_assigned_char(c::Char) Unicode.isassigned(c)
@deprecate normalize_string(s::AbstractString, nf::Symbol; kwargs...) Unicode.normalize(s, nf; kwargs...)
@deprecate normalize_string(s::AbstractString; kwargs...) Unicode.normalize(s; kwargs...)

# END 0.7 deprecations

end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI