text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# 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.
##########################################################################

# Dynamic Double Standard Model, with some modifications
# From paper: "A dynamic model and parallel tabu search heuristic for real-time ambulance relocation"

# Modifications (compared with original ddsm):
# - Added slack variables to coverage constraints, to avoid infeasibility. Cost of slack per unit is given by slackWeight.
# - Redeployment cost only depends on travel time, not on: ambulance history, avoiding round trips.

# initialise data relevant to move up
function initDdsm!(sim::Simulation;
	coverFractionTargetT1::Float = 0.5, travelTimeCost::Float = 50.0, slackWeight::Float = 1e9,
	coverTimeDemandPriorities::Vector{Priority} = [highPriority, lowPriority],
	options::Dict{Symbol,Any} = Dict{Symbol,Any}())
	# slackWeight is the weight to apply to slack variables, needs to be sufficiently large so that the slack variables are zero when possible
	# coverTimeDemandPriorities[i] is demand priority for coverTimes[i], where coverTimes[1] and [2] are the targets for ddsm

	@assert(0 <= coverFractionTargetT1 <= 1)
	@assert(length(coverTimeDemandPriorities) == 2)

	# initialise demand and demand coverage data if not already initialised
	sim.demand.initialised || initDemand!(sim)
	sim.demandCoverage.initialised || initDemandCoverage!(sim)
	demand = sim.demand # shorthand

	@assert(all([demandMode.rasterIndex for demandMode in demand.modes[demand.modeLookup[i,:]]] |> unique |> length == 1 for i = 1:demand.numSets)) # distribution of demand must be same for all demand priorities, in any given time period

	# two target cover times for ddsm
	coverTimes = [sim.demandCoverage.coverTimes[p] for p in coverTimeDemandPriorities]
	@assert(length(coverTimes) == 2)
	@assert(coverTimes[1] < coverTimes[2]) # change to warning?

	# set options
	ddsmd = sim.moveUpData.ddsmData # shorthand
	ddsmd.options[:solver] = "cbc" # can be slower than glpk, but more reliable for some reason
	ddsmd.options[:solver_args] = []
	ddsmd.options[:solver_kwargs] = []
	ddsmd.options[:v] = v"1"
	ddsmd.options[:z_var] = true
	ddsmd.options[:bin_tol] = 1e-5
	merge!(ddsmd.options, Dict([:x_bin => true, :y11_bin => true, :y12_bin => true, :y2_bin => true]))
	ddsmOptions!(sim, options)

	ddsmd.coverFractionTargetT1 = coverFractionTargetT1
	ddsmd.travelTimeCost = travelTimeCost
	ddsmd.slackWeight = slackWeight
	ddsmd.coverTimeDemandPriorities = coverTimeDemandPriorities # coverTimeDemandPriorities[i] is demand priority for coverTimes[i]
	ddsmd.coverTimes = coverTimes
end

# change the options for ddsm
function ddsmOptions!(sim, options::Dict{Symbol,T}) where T <: Any
	options = merge!(sim.moveUpData.ddsmData.options, options) # update options

	@assert(in(options[:solver], ["cbc", "glpk", "gurobi"]))
	@assert(typeof(options[:v]) == VersionNumber)
	@assert(typeof(options[:z_var]) == Bool)
	@assert(0 < options[:bin_tol] < 0.1)

	if options[:solver] == "gurobi" try Gurobi; catch; options[:solver] = "cbc"; @warn("Failed to use Gurobi, using Cbc instead.") end end

	# options[:v] == v"0" # do nothing, values should already be set in options if using this
	options[:v] == v"1" && merge!(options, Dict([:x_bin => true, :y11_bin => true, :y12_bin => true, :y2_bin => true]))
	options[:v] == v"2" && merge!(options, Dict([:x_bin => true, :y11_bin => true, :y12_bin => false, :y2_bin => false]))
	options[:v] == v"3" && merge!(options, Dict([:x_bin => false, :y11_bin => true, :y12_bin => false, :y2_bin => false]))
	options[:v] == v"4" && merge!(options, Dict([:x_bin => false, :y11_bin => true, :y12_bin => true, :y2_bin => true]))
	if options[:y11_bin] == false && in(options[:solver], ["cbc", "glpk"]) @warn("Removing binary constraint for `y11` may not work with CBC or GLPK.") end
	if options[:x_bin] == false && options[:solver] == "gurobi" @warn("Removing binary constraint for `x` may not work with Gurobi.") end
	@assert(all(key -> haskey(options, key), [:x_bin, :y11_bin, :y12_bin, :y2_bin]))

	if isdefined(sim, :backup) sim.backup.moveUpData.ddsmData.options = options end # to keep options if sim is reset

	return options
end

function ddsmMoveUp(sim::Simulation)
	@assert(sim.moveUpData.useMoveUp)

	# shorthand:
	@unpack ambulances, stations, numStations, demand = sim
	currentTime = sim.time
	ddsmd = sim.moveUpData.ddsmData
	@unpack coverTimes, coverTimeDemandPriorities, coverFractionTargetT1, slackWeight, options = ddsmd

	# get movable ambulances (movableAmbs)
	ambMovable = [isAmbMovable(amb) for amb in ambulances]
	movableAmbs = ambulances[ambMovable]
	numMovableAmbs = length(movableAmbs)

	if numMovableAmbs == 0
		return moveUpNull()
	end

	# calculate travel time for each movable ambulance
	ambToStationTimes = zeros(Float, numMovableAmbs, numStations) # ambToStationTimes[i,j] = time for movable ambulance i to travel to station j
	for i = 1:numMovableAmbs
		ambToStationTimes[i,:] = ambMoveUpTravelTimes!(sim, movableAmbs[i])
	end
	ambToStationCosts = ambToStationTimes * ddsmd.travelTimeCost

	# get demand point coverage data for the two target cover times
	# for coverTimes[ti], point j has demand pointDemands[ti][j] and is covered by stations pointStations[ti][j]
	pointStations = [] # pointStations[ti][j] is coverTimes[ti]
	pointDemands = []
	numPoints = Int[] # can have different number of points (really point sets) for each cover time
	demandPriorityArrivalRates = [getDemandMode!(demand, demandPriority, currentTime).arrivalRate for demandPriority in priorities]
	for demandPriority in coverTimeDemandPriorities
		# get demand point coverage data
		pointsCoverageMode = getPointsCoverageMode!(sim, demandPriority, currentTime)
		demandMode = getDemandMode!(demand, demandPriority, currentTime)
		pointSetsDemands = getPointSetsDemands!(sim, demandPriority, currentTime; pointsCoverageMode = pointsCoverageMode) * demandMode.rasterMultiplier
		pointSetsDemands *= sum(demandPriorityArrivalRates) / demandPriorityArrivalRates[Int(demandPriority)] # scale demand to be for all priorities, not just current demandPriority

		push!(pointStations, pointsCoverageMode.stationSets)
		push!(pointDemands, pointSetsDemands)
		push!(numPoints, length(pointSetsDemands))
	end
	@assert(isapprox(sum(pointDemands[1]), sum(pointDemands[2])))
	# could have it so that pointStations, pointDemands, and numPoints are only updated if the demandSet has changed since last query,
	# but this may only shave off a couple of seconds per 100 sim days, which is not the current bottleneck

	######################
	# IP

	# shorthand
	a = numMovableAmbs
	np = numPoints # np[ti] is number of points for coverTimes[i]
	s = numStations
	t = coverTimes

	# using JuMP
	model = Model()
	solver = options[:solver] # shorthand
	args = options[:solver_args] # shorthand
	kwargs = options[:solver_kwargs] # shorthand
	if solver == "cbc" set_optimizer(model, with_optimizer(Cbc.Optimizer, logLevel=0, args...; kwargs...))
	elseif solver == "glpk" set_optimizer(model, with_optimizer(GLPK.Optimizer, args...; kwargs...))
	elseif solver == "gurobi" @stdout_silent(set_optimizer(model, with_optimizer(Gurobi.Optimizer, OutputFlag=0, args...; kwargs...)))
	end

	m = model # shorthand
	options[:x_bin] ? @variable(m, x[i=1:a,j=1:s], Bin) : @variable(m, 0 <= x[i=1:a,j=1:s] <= 1)
	options[:y11_bin] ? @variable(m, y11[p=1:np[1]], Bin) : @variable(m, 0 <= y11[p=1:np[1]] <= 1)
	options[:y12_bin] ? @variable(m, y12[p=1:np[1]], Bin) : @variable(m, 0 <= y12[p=1:np[1]] <= 1)
	options[:y2_bin] ? @variable(m, y2[p=1:np[2]], Bin) : @variable(m, 0 <= y2[p=1:np[2]] <= 1)
	# x[i,j] = 1 if ambulance: movableAmbs[i] should be moved to station: stations[j], 0 otherwise
	# y11[p,k] = 1 if demand point p is covered at least once within t[1], 0 otherwise
	# y12[p,k] = 1 if demand point p is covered at least twice within t[1], 0 otherwise
	# y2[p] = 1 if demand point p is covered at least once within t[2], 0 otherwise

	@variables(model, begin
		(s1 >= 0) # slack
		(s2 >= 0) # slack
	end)

	@constraints(model, begin
		(ambAtOneStation[i=1:a], sum(x[i,:]) == 1) # each ambulance must be assigned to one station
		(pointCoverOrderY1[p=1:np[1]], y11[p] >= y12[p]) # single coverage before double coverage; not needed for y2
		(demandCoveredOnceT1, sum(y11[p] * pointDemands[1][p] for p=1:np[1]) + s1 >= coverFractionTargetT1 * sum(pointDemands[1])) # fraction of demand covered once within t[1]
		(demandCoveredOnceT2, sum(y2[p] * pointDemands[2][p] for p=1:np[2]) + s2 >= sum(pointDemands[2])) # all demand covered once within t[2]
	end)

	if options[:z_var]
		@variable(model, z[j=1:s], Int) # z[j] = number of ambulances to assign to station j
		@constraints(model, begin
			(stationAmbCount[j=1:s], z[j] == sum(x[:,j]))
			(pointCoverCountY1[p=1:np[1]], y11[p] + y12[p] <= sum(z[pointStations[1][p]]))
			(pointCoverCountY2[p=1:np[2]], y2[p] <= sum(z[pointStations[2][p]]))
		end)
	else
		@constraints(model, begin
			(pointCoverCountY1[p=1:np[1]], y11[p] + y12[p] <= sum(x[:,pointStations[1][p]]))
			(pointCoverCountY2[p=1:np[2]], y2[p] <= sum(x[:,pointStations[2][p]]))
		end)
	end

	@expressions(model, begin
		demandCoveredTwiceT1, sum(y12[p] * pointDemands[1][p] for p=1:np[1])
		totalAmbTravelCost, sum(x[i,j] * ambToStationCosts[i,j] for i=1:a, j=1:s)
		slackCost, (s1 + s2) * slackWeight
	end)

	# solve
	@objective(model, Max, demandCoveredTwiceT1 - totalAmbTravelCost - slackCost)
	@stdout_silent optimize!(model)
	@assert(termination_status(model) == MOI.OPTIMAL)

	# get solution
	vals = Dict()
	vals[:x] = JuMP.value.(x)
	vals[:y11] = JuMP.value.(y11)
	vals[:y12] = JuMP.value.(y12)
	vals[:y2] = JuMP.value.(y2)

	sol = convert(Array{Bool,2}, round.(vals[:x]))

	if checkMode
		@assert(all(sum(sol, dims=2) .== 1)) # check constraint: ambAtOneStation
		# check that values are binary
		for sym in [:x, :y11, :y12, :y2]
			err = maximum(abs.(vals[sym] - round.(vals[sym])))
			@assert(err <= options[:bin_tol], (sym, err))
		end
	end

	ambStations = [stations[findfirst(sol[i,:])] for i=1:a]

	return movableAmbs, ambStations
end


text:
#julia
using Makie
scene = Scene()
N = 32
function xy_data(x,y,i, N)
    x = ((x/N)-0.5f0)*i
    y = ((y/N)-0.5f0)*i
    r = sqrt(x*x + y*y)
    res = Float32(sin(r)/r)
    isnan(res) ? 1f0 : res
end

surf_func(i) = [Float32(xy_data(x, y, i, 32)) + 0.5 for x=1:32, y=1:32]

z = surf_func(20)
range = linspace(0, 3, N)
surf = surface(range, range, z, colormap = :Spectral)
wf = wireframe(range, range, surf[:z] .+ 1.0,
    linewidth = 2f0, color = lift_node(x-> x[5], surf[:colormap])
)
axis(linspace(0, 3, 4), linspace(0, 3, 4), linspace(0, 3, 4))
center!(scene)

wf[:linewidth] = 1
surf[:colormap] = :YlGnBu
io = VideoStream(scene)
for i in linspace(0, 60, 200)
    surf[:z] = surf_func(i)
    recordframe!(io)
end
io

ext: jl


text:
# Compare performance of a duck-typed function with trait-dispatched function

## a simple function:
@traitfn function ft2{X,Y; Arith{X,Y}}(x::X,y::Y)
    out = zero(promote(x,y)[1])
    for xe in 1:round(Int,x)
        out += xe + y
    end
    out
end
function ff2{X,Y}(x::X,y::Y)
    out = zero(promote(x,y)[1])
    for xe in 1:round(Int,x)
        out += xe + y
    end
    out
end

n = 100000
x = rand(n)
y = rand(n)

function gt()
    out = 0.0
    for i=1:n
        ft2(x[i],y[i])
    end
    out
end
function gg()
    out = 0.0
    for i=1:n
        ff2(x[i],y[i])
    end
    out
end

gg()
gc()
t, m = @timeit gg()
println("Duck-typed function:  time $(round(t,5))s, mem $(round(Int,m/1e6))MB")

gt()
gc()
t, m = @timeit gt()
println("Trait-typed function: time $(round(t,5))s, mem $(round(Int,m/1e6))MB")
##


text:
function floydwarshall_test()
    A = load_matrix_network("all_shortest_paths_example")
    #A.nzval += abs(minimum(A)) + 1 # remove negative edges
    nzvals = nonzeros(A)
    val = abs(minimum(A)) + 1
    for i=1:length(nzvals)
        nzvals[i] += val
    end
    m = size(A,1)
    D2 = zeros(Float64,m,m)
    P2 = zeros(Float64,m,m)
    for i=1:m
        (d,p) = dijkstra(A,i)
        D2[i,:] = d
        P2[i,:] = p
    end
    (D,P) = floydwarshall(A)

    if !isequal(D,D2)
        error("Floyd Warshall failed: Incorrect distances")
    end

    if !isequal(P,P2)
        error("Floyd Warshall failed: Incorrect predecessors")
    end
    return true
end


text:
module RLEnvironments
    using Distributions: Uniform
    import Distributions.sample

    abstract type Environment{S,A,R<:Real} end
    abstract type DiscEnv{S,A,R} <: Environment{S,A,R} end
    abstract type ContEnv{S,A,R} <: Environment{S,A,R} end
    abstract type DiffEnv{S,A,R} <: Environment{S,A,R} end
    function step! end
    function reset! end

    #DiscEnv: Vector{A}, ContEnv: Vector{Tuple{A, A}}
    struct ActionSpace{A}
        actions::Vector{A}
    end

    ActionSpace(n::A) where A <: Real = ActionSpace(collect(one(A):n))

    Base.length(a::ActionSpace) = length(a.actions)
    Base.eltype(::Type{ActionSpace{A}}) where A = A
    Base.eltype(::Type{ActionSpace{Tuple{A, A}}}) where A = A

    sample(a::ActionSpace) = rand(a.actions)
    function sample(a::ActionSpace{Tuple{T, T}}) where T <: AbstractFloat
        return map(x -> T(rand(Uniform(x[1], x[2]))), a.actions)
    end

    #n = numstates // n = length(state)
    struct ObservationSpace{S}
        n::Int64
    end

    Base.length(o::ObservationSpace) = o.n
    Base.eltype(::Type{ObservationSpace{S}}) where S = S

    function test(f, env::Environment{S, A, R}, maxt=100) where {S, A, R}
        reward = zero(R)
        s = reset!(env)
        for _ in 1:maxt
            s, r, done = step!(env, f(s))
            reward += r
            done && break
        end
        return reward
    end

    include("bandit.jl")
    include("simplerooms.jl")
    include("mountaincar.jl")
    include("cartpole.jl")
    include("pendulum.jl")

    export Environment, DiscEnv, ContEnv, DiffEnv, step!, reset!, test
    export Bandit, SimpleRooms, CartPole, CartPoleContinuous, Pendulum,
        MountainCar, MountainCarContinuous
end


text:
using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")
function sigmoid_(z)

    return 1.0 ./ (1.0 .+ exp.(-z))

end

H = x->sigmoid_((x-1)*100)

x = LinRange(0,2,1000)
y = H.(x)
plot(x, y, linewidth=4)
plot(ones(100), LinRange(0,1,100), "--", linewidth=4)
text(1.1, 0.1, "\$\\tilde\\sigma_Y\$", fontsize=20)

text(1.1, 0.6, "\$\\sigma_{\\mathrm{plasticity}}\$",  fontsize=20)
text(1.1, 0.5, "dominates",  fontsize=20)

text(0.1, 0.6, "\$\\sigma_{\\mathrm{elasticity}}\$", fontsize=20)
text(1.1, 0.5, "dominates",  fontsize=20)

xlabel("\$\\sigma\$")
ylabel("\$D(\\sigma, \\tilde\\sigma_Y)\$")
mpl.save("H.tex")


text:
module BNPRegressionGGA2021

using BayesNegativeBinomial
using Distributions
using HypergeometricFunctions
using LinearAlgebra
using Random
using SpecialFunctions
using StatsBase: mean_and_std
using Walker2020Sampling
import Distributions: pdf, logpdf

const BNB = BayesNegativeBinomial

# All DGPM models are subtypes of AbstractModel:
abstract type AbstractModel end

# All subtypes of AbstractModel must extend (by composition) the Skeleton type:
include("skeleton.jl")

# All subtypes of AbstractModel must implement the following interface:
include("interface.jl")

# The methods available for any AbstractModel are described here:
include("methods.jl")

# The specific models are implemented here:
include("dgsbpnormaldependent.jl")
include("dgsbpnormal.jl")
include("dgpmerlang.jl")
include("dgpmbeta.jl")

end


text:
function tgenerator(dt,nt,ns,fold)
j=1;
J=zeros(ns);
for k=1:ns
    J[k]=j;
    j=j+floor(nt/fold+((rand(1)[1])*(nt/fold)));
end
tt=dt.*J';
end

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
export  execute,
        Environment,
        Action, Action_NoOp,
        Percept,
        Agent,
        AgentProgram,
            TableDrivenAgentProgram,
            ReflexVacuumAgentProgram,
            SimpleReflexAgentProgram,
            ModelBasedReflexAgentProgram,
        State,
        Rule

"""
An agent perceives an *environment* through sensors and acts with actuators.

Sensors provide agent the *percepts*, based on which the agent delivers
*actions*

Pg. 35, AIMA 3ed
"""
abstract type Environment end

"""
*AgentProgram* is an internal representation of an agent function with a
concrete implementation. While *agent function* can be abstract *AgentProgram*
provides clear direction to the implementation.

Pg. 35, AIMA 3ed
"""
abstract type AgentProgram end

"""
*Action* is an agent's response to the environment through actuators.

Although the representation of a string may suffice for more most sample
programs, an abstract type is introduced to emphasize the need for
providing a concrete type based on the environment or agent at hand.

In most problems we try to solve, the *Action* may be driven by the choice
of the *Environment*
"""
abstract type Action end

"""
*NoOp* is a directive where the agent does not take any further action.
"""
struct NoOpActionType <: Action
  val::Symbol
  NoOpActionType() = new(Symbol("NoOp"))
end

const Action_NoOp = NoOpActionType()

"""
*Percept* is an input to the *Agent* from environment through sensors.

Although the representation of a Tuple may suffice for more most sample
programs, an abstract type is introduced to emphasize the need for
providing a concrete type based on the environment or agent at hand.

In most problems we try to solve, the *Percept* may be driven by the choice
of the *Environment*
"""
abstract type Percept end

"""
Given a *Percept* returns an *Action* apt for the agent.

Depending on the agent program the function may respond with different *Action*
evaluation strategies.
"""
execute(ap::AP, p::Percept) where {AP <: AgentProgram} = error(E_ABSTRACT)

"""
*TableDrivenAgentProgram* is a simple model of an agent program where all
percept sequences are well-known ahead in time and can be organized as a
mapping from percepts to action.

Look at the corresponding execute method for *Action* evaluation strategy.

The implementation must have the following methods:

1. append - percept to the list of percepts seen my the AgentProgram
2. lookup - the percepts in the tables of the AgentProgram

Fig 2.7 Pg. 47, AIMA 3ed
"""
abstract type TableDrivenAgentProgram <: AgentProgram end

function execute(ap::TableDrivenAgentProgram, percept::Percept)
    append(ap.percepts, percept)
    action = lookup(ap.table, ap.percepts)
    return action
end

"""
*Rule* is an abstract representation of a framework that associates a *State*
condition to the appropriate action.

Definition of a condition can be implementation dependent.
"""
abstract type Rule end

"""
*State* is an internal evaluated position of the Environment. In the context
of the problem the *Environment* can be one of the stated states. Any input or'
action may lead to change in *Environment* state.
"""
abstract type State end

"""
*SimpleReflexAgentProgram* is a simple *Percept* to *Action* matching state
based rules.

It does not depend on the historical percept data.

It needs to implement two methods

1. interpret_input
2. rule_match

for all the concrete implementations.

3. rules - Will provide all the rules associated with the
AgentProgram.

Fig 2.10 Pg. 49, AIMA 3ed
"""
abstract type SimpleReflexAgentProgram <: AgentProgram end

function execute(ap::SimpleReflexAgentProgram, percept::Percept)
    state = interpret_input(percept);
    rule = rule_match(state, ap.rules);
    action = rule.action;
    return action;
end

interpret_input(percept) = error(E_ABSTRACT)
rule_match(state, rules) = error(E_ABSTRACT)

"""
*ModelBasedReflexAgentProgram* uses a model which is close to the
understanding of the world.

The *AgentProgram* updates the states based on the *Percepts* received.
"""
abstract type ModelBasedReflexAgentProgram <: AgentProgram end

function execute(ap::ModelBasedReflexAgentProgram, percept::Percept)
    ap.state = update_state(ap.state, ap.action, percept, ap.model)
    rule = rule_match(ap.state, ap.rules)
    ap.action = rule.action
    return ap.action
end

update_state(state, action, percept, model) = error(E_ABSTRACT)

"""
*SimpleProblemSolvingAgentProgram* formulates the goal then the problem.
Searches for the sequence of *Actions* that would solve the problem.

pg.67 AIMA 3e
"""
abstract type SimpleProblemSolvingAgentProgram <: AgentProgram end

function execute(ap::SimpleProblemSolvingAgentProgram, percept::Percept)
    #persistent: seq, an action sequence, initially empty
    #state, some description of the current world state
    #goal , a goal, initially null
    #problem, a problem formulation
    ap.state = update_state(ap.state, percept)
    if isempty(ap.seq)
        ap.goal = formulate_goal(ap.state)
        ap.problem = formulate_problem(ap.state, ap.goal)
        ap.seq = search(ap.problem)
        if ap.seq == failure
            return nothing
        end
    end
    action = first(ap.seq)
    ap.seq = rest(ap.seq)
    return action
end

update_state(state, percept) = error(E_ABSTRACT)
formulate_goal(state) = error(E_ABSTRACT)
formulate_problem(state, goal) = error(E_ABSTRACT)

"""
Agent perceives *Environment* through sensors and acts based on actuators.

While this code may not be there in the book this is a general outline one can
gather from the description in the book.

pg. 35 Fig.2.1 AIMA 3e
"""
struct Agent{AP}
    program::AP
    Agent{AP}(ap::AP) where{AP <: AgentProgram} = new(ap)
end

execute(a::Agent, percept::Percept) = execute(a.program, percept)


text:
# Compute the gradient of 2D array, returns grad = [real, imag]

function gradient_2D(arr, X)
    arr_real = real(arr)
    arr_imag = imag(arr)
    x = X[1]; y = X[2];
    dx = x[2] - x[1];
    dy = y[2] - y[1];
    dxarr_real = similar(arr_real);
    dyarr_real = similar(arr_real);
    dxarr_imag = similar(arr_imag);
    dyarr_imag = similar(arr_imag);

    for xin in 2:length(x)
        dxarr_real[xin, :] .= (arr_real[xin, :] .- arr_real[xin-1, :]) ./ dx;
        dxarr_imag[xin, :] .= (arr_imag[xin, :] .- arr_imag[xin-1, :]) ./ dx;
    end
    dxarr_real[1, :] .= 0;
    dxarr_imag[1, :] .= 0;

    for yin in 2:length(y)
        dyarr_real[:, yin] .= (arr_real[:, yin] .- arr_real[:, yin-1]) ./ dy;
        dyarr_imag[:, yin] .= (arr_imag[:, yin] .- arr_imag[:, yin-1]) ./ dy;
    end
    dyarr_real[:, 1] .= 0;
    dyarr_imag[:, 1] .= 0;
    GRADarr_real = [dxarr_real, dyarr_real];
    GRADarr_imag = [dxarr_imag, dyarr_imag];
    grad = [GRADarr_real, GRADarr_imag];
    return grad;
end

# Compute the gradient of 3D array, returns grad = [real, imag]
function gradient_3D(arr, X)

    arr_real = real(arr)
    arr_imag = imag(arr)
    x = X[1]; y = X[2]; z = X[3];
    dx = x[2] - x[1];
    dy = y[2] - y[1];
    dz = z[2] - z[1];

    dxarr_real = similar(arr_real);
    dyarr_real = similar(arr_real);
    dzarr_real = similar(arr_real);
    dxarr_imag = similar(arr_imag);
    dyarr_imag = similar(arr_imag);
    dzarr_imag = similar(arr_imag);

    for xin in 2:length(x)
        dxarr_real[xin, :, :] .= (arr_real[xin, :, :] .- arr_real[xin-1, :, :]) ./ dx;
        dxarr_imag[xin, :, :] .= (arr_imag[xin, :, :] .- arr_imag[xin-1, :, :]) ./ dx;
    end
    dxarr_real[1, :, :] .= 0;
    dxarr_imag[1, :, :] .= 0;

    for yin in 2:length(y)
        dyarr_real[:, yin, :] .= (arr_real[:, yin, :] .- arr_real[:, yin-1, :]) ./ dy;
        dyarr_imag[:, yin, :] .= (arr_imag[:, yin, :] .- arr_imag[:, yin-1, :]) ./ dy;
    end
    dyarr_real[:, 1, :] .= 0;
    dyarr_imag[:, 1, :] .= 0;

    for zin in 2:length(z)
        dzarr_real[:, :, zin] .= (arr_real[:, :, zin] .- arr_real[:, :, zin - 1]) ./ dz;
        dzarr_imag[:, :, zin] .= (arr_imag[:, :, zin] .- arr_imag[:, :, zin - 1]) ./ dz;
    end
    dzarr_real[:, :, 1] .= 0;
    dzarr_imag[:, :, 1] .= 0;

    GRADarr_real = [dxarr_real, dyarr_real, dzarr_real];
    GRADarr_imag = [dxarr_imag, dyarr_imag, dzarr_imag];
    GRADpsi = [GRADarr_real, GRADarr_imag];
    return GRADpsi;
end

# Pseudo-vorticity 3D, needs X = [x, y, z] for the grid and array gradient as parameters.
# Returns pseudo-vorticity vector
function wps_v(gradpsi, X)
    x = X[1]; y = X[2]; z = X[3];
    grad_real = gradpsi[1]; grad_imag = gradpsi[2];
    dx_r = grad_real[1]; dy_r = grad_real[2]; dz_r = grad_real[3];
    dx_i = grad_imag[1]; dy_i = grad_imag[2]; dz_i = grad_imag[3];
    grad_v_r = [dx_r[x,y,z], dy_r[x,y,z], dz_r[x,y,z]];
    grad_v_i = [dx_i[x,y,z], dy_i[x,y,z], dz_i[x,y,z]];

    cross_v = cross(grad_v_r,  grad_v_i)
    wps = cross_v
    return wps;
end

# Pseudo-vorticity 2D, needs X = [x, y] for the grid and array gradient as parameters.
# Returns vector orthog to x, y. Good to sanity check as wps = 0 if no nearby vorticies.
function wps_v_2D(gradpsi, X)
    x = X[1]; y = X[2];
    grad_real = gradpsi[1]; grad_imag = gradpsi[2];
    dx_r = grad_real[1]; dy_r = grad_real[2];
    dx_i = grad_imag[1]; dy_i = grad_imag[2];
    grad_v_r = [dx_r[x,y], dy_r[x,y], 0];
    grad_v_i = [dx_i[x,y], dy_i[x,y], 0];

    cross_v = cross(grad_v_r, grad_v_i)
    wps = cross_v;
    return wps;
end

# Find closest element to x in array
function find_nearest(x, arr)
    nearest = arr[1];
    nearest_idx = 1;
    diff = abs(x-nearest);

    for i in 2:length(arr)
        if abs(x-arr[i]) < diff
            nearest = arr[i];
            nearest_idx = i;
            diff = abs(x-arr[i]);
        end
    end
    return [nearest,  Int(nearest_idx)];
end

function findidx_uniform(x, arr)
    # For arrays where xi is between [-L, L] and uniform spacing
    # So xi = -L/2 + (L/N)*i
    # => i = (N/L)*(x + L/2)
    N = length(arr);
    L = arr[end]*2;
    idx = (N/L)*(x + L/2);
    if idx < 1
        println("x is not within the bounds of this array");
        return;
    else
        return idx;
    end
end

function wps_array_plot(psi, zslice)
    grad = gradient_3D_cent(psi, X);
    psi_torus = Torus(psi[:, :, zslice], x, y);
    vorts = vortex_array(findvortices(psi_torus))
    wps_arr = similar(real(psi[:, :, 1]))

    for i in 1:length(x)
        for j in 1:length(y)
            wps_arr[i, j] = norm(wps_v(grad, [i,j, zslice]))
        end
    end
    wps_arr = abs.(wps_arr)
    max = findmax(wps_arr)[1]
    wps_arr = wps_arr ./max

    heatmap(x, y, wps_arr)
    scatter!(vorts[:, 2], vorts[:, 1])
end

function gradient_3D_cent(arr, X)
    x = X[1]; y = X[2]; z = X[3];
    arr_r = real(arr);
    arr_i = imag(arr);

    dx = x[2] - x[1]; dy = y[2] - y[1]; dz = z[2] - z[1];
    dxarr_r = similar(arr_r);
    dyarr_r = similar(arr_r);
    dzarr_r = similar(arr_r);
    dxarr_i = similar(arr_i);
    dyarr_i = similar(arr_i);
    dzarr_i = similar(arr_i);

    # Periodic Boundary conditions
    dxarr_r[1, :, :] .= (arr_r[2, :, :] .- arr_r[end, :, :]) ./ (2*dx);
    dxarr_r[end, :, :] .= (arr_r[1, :, :] .- arr_r[end-1, :, :]) ./ (2*dx);
    dxarr_i[1, :, :] .= (arr_i[2, :, :] .- arr_i[end, :, :]) ./ (2*dx);
    dxarr_i[end, :, :] .= (arr_i[1, :, :] .- arr_i[end-1, :, :]) ./ (2*dx);

    for i in 2:length(x)-1
        dxarr_r[i,:,:] .= (arr_r[i+1, :, :] .- arr_r[i-1, :, :]) ./ (2*dx);
        dxarr_i[i,:,:] .= (arr_i[i+1, :, :] .- arr_i[i-1, :, :]) ./ (2*dx);
    end

    dyarr_r[:, 1, :] .= (arr_r[:, 2, :] .- arr_r[:, end, :]) ./ (2*dy);
    dyarr_r[:, end, :] .= (arr_r[:, 1, :] .- arr_r[:, end-1, :]) ./ (2*dy);
    dyarr_i[:, 1, :] .= (arr_i[:, 2, :] .- arr_i[:, end, :]) ./ (2*dy);
    dyarr_i[:, end, :] .= (arr_i[:, 1, :] .- arr_i[:, end-1, :]) ./ (2*dy);

    for i in 2:length(y)-1
        dyarr_r[:,i,:] .= (arr_r[:, i+1, :] .- arr_r[:, i-1, :]) ./ (2*dy);
        dyarr_i[:,i,:] .= (arr_i[:, i+1, :] .- arr_i[:, i-1, :]) ./ (2*dy);
    end

    dzarr_r[:, :, 1] .= (arr_r[:, :, 2] .- arr_r[:, :, end])./ (2*dz);
    dzarr_r[:, :, end] .= (arr_r[:, :, 1] .- arr_r[:, :, end-1])./ (2*dz);
    dzarr_i[:, :, 1] .= (arr_i[:, :, 2] .- arr_i[:, :, end])./ (2*dz);
    dzarr_i[:, :, end] .= (arr_i[:, :, 1] .- arr_i[:, :, end-1])./ (2*dz);

    for i in 2:length(z)-1
        dzarr_r[:, :, i] .= (arr_r[:, :, i+1] - arr_r[:, :, i-1]) ./ (dz*2);
        dzarr_i[:, :, i] .= (arr_i[:, :, i+1] - arr_i[:, :, i-1]) ./ (dz*2);
    end

    grad_r = [dxarr_r, dyarr_r, dzarr_r];
    grad_i = [dxarr_i, dyarr_i, dzarr_i];
    grad = [grad_r, grad_i];

    return grad;
end

function grad_itp(grad, X)
    x = X[1]; y = X[2]; z = X[3];
    x = LinRange(x[1], x[end], length(x));
    y = LinRange(y[1], y[end], length(y));
    z = LinRange(z[1], z[end], length(z));
    grad_r = grad[1]; grad_i = grad[2];
    grx = grad_r[1];
    gry = grad_r[2];
    grz = grad_r[3];
    gix = grad_i[1];
    giy = grad_i[2];
    giz = grad_i[3];
    grx_itp = interpolate(grx, BSpline(Cubic(Line(OnGrid()))));
    gry_itp = interpolate(gry, BSpline(Cubic(Line(OnGrid()))));
    grz_itp = interpolate(grz, BSpline(Cubic(Line(OnGrid()))));
    gix_itp = interpolate(gix, BSpline(Cubic(Line(OnGrid()))));
    giy_itp = interpolate(giy, BSpline(Cubic(Line(OnGrid()))));
    giz_itp = interpolate(giz, BSpline(Cubic(Line(OnGrid()))));
    sgrx_itp = scale(grx_itp, x, y, z);
    sgry_itp = scale(gry_itp, x, y, z);
    sgrz_itp = scale(grz_itp, x, y, z);
    sgix_itp = scale(gix_itp, x, y, z);
    sgiy_itp = scale(giy_itp, x, y, z);
    sgiz_itp = scale(giz_itp, x, y, z);
    gr_v = [sgrx_itp, sgry_itp, sgrz_itp];
    gi_v = [sgix_itp, sgiy_itp, sgiz_itp];
    grad = [gr_v, gi_v];
    return grad;
end

function wps_int(grad, Xin)

    grad_r = grad[1]; grad_i = grad[2];

    gr_v = [grad_r[1](Xin[1], Xin[2], Xin[3]), grad_r[2](Xin[1], Xin[2], Xin[3]), grad_r[3](Xin[1], Xin[2], Xin[3])];
    gi_v = [grad_i[1](Xin[1], Xin[2], Xin[3]), grad_i[2](Xin[1], Xin[2], Xin[3]),  grad_i[3](Xin[1], Xin[2], Xin[3])];

    wps = cross(gr_v, gi_v);
    return wps;
end

function psi_itp(psi, X)
    x = X[1]; y = X[2]; z = X[3];
    x = LinRange(x[1], x[end], length(x));
    y = LinRange(y[1], y[end], length(y));
    z = LinRange(z[1], z[end], length(z));
    psi_itp = interpolate(psi, BSpline(Cubic(Line(OnGrid()))));
    psi_itp = scale(psi_itp, x, y, z);

    return psi_itp
end

function vortices_plot(psi, x, y, zidx)
    psi = Torus(psi_xspc[:, :, zidx], x, y);
    vort_array = vortex_array(findvortices(psi));
    vort_pos = vort_array[:, 1:2];

    @views Plots.heatmap(x, y, angle.(psi_xspc[:, :, zidx]), title = "z level = $zidx", label = "Phase")
    Plots.scatter!(vort_pos[:, 2], vort_pos[:, 1], color =:white, label = "Vortices")
    Plots.xlabel!("x")
    Plots.ylabel!("y")

    #when calling:
    #=
    gr()
    anim = @animate for zidx = 1:Nz
        vortices_plot(psi_kspc, x, y, zidx)
    end
    gif(anim, "anim_fps15.gif", fps = 8)
    =#

end

function find_closest_tuple(arr, X)
    x = X[1]; y = X[2];
    closest_dist = 1000;
    index = 0;
    for i in 1:length(arr[:, 1])
        dist = (x - arr[i, 1])^2 + (y - arr[i, 2])^2;
        dist = sqrt(dist);
        if dist < closest_dist
            index = i;
            closest_dist = dist;
        end
    end

    return index;
end


text:
#=

   From https://math.stackexchange.com/questions/513250/conditional-probability-question-with-cards-where-the-other-side-color-is-to-be
   """
   A box contains three cards. One card is red on both sides, one card is green on both sides,
   and one card is red on one side and green on the other. One card is selected from the
   box at random, and the color on one side is observed. If this side is green, what is
   the probability that the other side of the card is also green?

   ... the answer to this question is 2/3.
   """

   Give that the color of card is green, then the probability
   of the different cards are:

   Distributions of variable selected_card
   all_green  =>   67030  (0.670300)
   red_green  =>   32970  (0.329700)

   I.e. the probability that it's the all green card is 2/3.

   Cf ~/webppl/card_problem.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function card_problem()
    red = 1
    green = 2

    all_red = 1
    all_green = 2
    red_green = 3

    #  What card did we select?
    selected_card ~ Categorical(simplex([1/3,1/3,1/3])) # [all_red,all_green,red_green]

    #  What is the color of the card (one side) that we see?
    card ~
        selected_card == all_red   ?  Dirac(red) :
        selected_card == all_green ?  Dirac(green) :
        # [red,green]
        selected_card == red_green ?  Categorical([1/2,1/2]) : Dirac(0)

    # Alternative:
    # if selected_card == all_red
    #     card ~  Dirac(red)
    # elseif selected_card == all_green
    #     card ~ Dirac(green)
    # else
    #     card ~ Categorical([1/2,1/2])
    # end

    #  The color of the card we see is green
    true ~ Dirac(card == green)

end

model = card_problem()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(10_000), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,10), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :selected_card,["all_red","all_green","red_green"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using PEG
@rule grammar = "using PEG\n" & rule[*]
@rule rule = r"@rule"p & nonterminal & r"="p & choice
@rule choice = seq & (r","p & seq)[*]
@rule seq = item & (r"&"p & item)[*] & (r"\|?>"p & julia_function)[:?]
@rule item = lookahead , counted
@rule lookahead = r"\("p & (r"[+-]"p) & seq & r"\)"p
@rule counted = single & (count)[:?]
@rule count = range , r"\["p & (":?" , r"[\*\+]"p) & r"]"p
@rule range = r"\["p & integer & (r":"p & (integer , r"end"w))[:?] & r"]"p
@rule integer = r"\d+"w
@rule single = parens , terminal , nonterminal
@rule parens = r"\("p & choice & r"\)"p
@rule nonterminal = r"\pL\w+"w
@rule terminal = regex , string & r"\s*"
@rule regex = r"\br" & string & r"[himpswx]*\s*"
@rule string = r"\"(\\.|[^\"])*\""
@rule julia_function = nonterminal
using Test
@test grammar(replace(open(s->read(s, String), @__FILE__), r"using Test.*"s => ""))[2] == ""


text:
# Learning with knowledge

export guess_example_value, generate_powerset, current_best_learning,
        version_space_learning,
        is_consistent_determination, minimal_consistent_determination,
        FOILKnowledgeBase, extend_example, choose_literal, new_literals,
        new_clause, foil;

function disjunction_value(e::Dict, d::Dict)
    for (k, v) in d
        if (!(typeof(v) <: AbstractString))
            error("disjunction_value(): Found an unexpected type, ", typeof(v), "!");
        end
        # Check for negation
        if (v[1] == '!')
            if (e[k] == v[2:end])
                return false;
            end
        elseif (e[k] != v)
            return false;
        end
    end
    return true;
end

"""
    guess_example_value(e::Dict, h::AbstractVector)

Return a guess for the logical value of the given example 'e' based on the given hypothesis 'h'.
"""
function guess_example_value(e::Dict, h::AbstractVector)
    for d in h
        if (disjunction_value(e, d))
            return true;
        end
    end
    return false;
end

function example_is_consistent(e::Dict, h::AbstractVector)
    return (e["GOAL"] == guess_example_value(e, h));
end

function example_is_false_positive(e::Dict, h::AbstractVector)
    if (e["GOAL"] == false)
        if (guess_example_value(e, h))
            return true;
        end
    end
    return false;
end

function example_is_false_negative(e::Dict, h::AbstractVector)
    if (e["GOAL"] == true)
        if (!(guess_example_value(e, h)))
            return true;
        end
    end
    return false;
end

function check_all_consistency(examples::AbstractVector, h::AbstractVector)
    for example in examples
        if (!(example_is_consistent(example, h)))
            return false;
        end
    end
    return true;
end

function specializations(prior_examples::AbstractVector, h::AbstractVector)
    local hypotheses::AbstractVector = [];
    for (i, disjunction) in enumerate(h)
        for example in prior_examples
            for (k, v) in example
                if ((haskey(disjunction, k)) || k == "GOAL")
                    continue;
                end

                local h_prime::Dict = copy(h[i]);
                h_prime[k] = "!" * v;
                local h_prime_prime::AbstractVector = copy(h);
                h_prime_prime[i] = h_prime;

                if (check_all_consistency(prior_examples, h_prime_prime))
                    push!(hypotheses, h_prime_prime);
                end
            end
        end
    end
    shuffle!(RandomDeviceInstance, hypotheses);
    return hypotheses;
end

function check_negative_consistency(examples::AbstractVector, h::Dict)
    for example in examples
        if (example["GOAL"])
            continue;
        end
        if (!example_is_consistent(example, [h]))
            return false;
        end
    end
    return true;
end

function generate_powerset(array::AbstractVector)
    local result::AbstractVector = Array{Any, 1}([()]);
    for element in array
        for i in eachindex(result)
            push!(result, (result[i]..., element));
        end
    end
    return Set{Tuple}(result);
end

function add_or_examples(prior_examples::AbstractVector, h::AbstractVector)
    local result::AbstractVector = [];
    local example::Dict = prior_examples[end];
    local attributes::Dict = Dict((k, v) for (k, v) in example if (k != "GOAL"));
    local attribute_powerset = setdiff!(generate_powerset(collect(keys(attributes))), Set([()]));

    for subset in attribute_powerset
        local h_prime::Dict = Dict();
        for key in subset
            h_prime[key] = attributes[key];
        end
        if (check_negative_consistency(prior_examples, h_prime))
            local h_prime_prime::AbstractVector = copy(h);
            push!(h_prime_prime, h_prime);
            push!(result, h_prime_prime);
        end
    end

    return result;
end

function generalizations(prior_examples::AbstractVector, h::AbstractVector)
    local hypotheses::AbstractVector = [];
    # Remove the empty set from the powerset.
    local disjunctions_powerset::Set = setdiff!(generate_powerset(collect(1:length(h))), Set([()]));
    for disjunctions in disjunctions_powerset
        h_prime = copy(h);
        deleteat!(h_prime, disjunctions);

        if (check_all_consistency(prior_examples, h_prime))
            append!(hypotheses, h_prime);
        end
    end

    for (i, disjunction) in enumerate(h)
        local attribute_powerset::Set = setdiff!(generate_powerset(collect(keys(disjunction))), Set([()]));
        for attributes in attribute_powerset
            h_prime = copy(h[i]);

            if (check_all_consistency(prior_examples, [h_prime]))
                local h_prime_prime::AbstractVector = copy(h);
                h_prime_prime[i] = copy(h_prime);
                push!(hypotheses, h_prime_prime);
            end
        end
    end
    if ((length(hypotheses) == 0) || (hypotheses == [Dict()]))
        hypotheses = add_or_examples(prior_examples, h);
    else
        append!(hypotheses, add_or_examples(prior_examples, h));
    end

    shuffle!(hypotheses);
    return hypotheses;
end

"""
    current_best_learning(examples::AbstractVector, h::AbstractVector, prior_examples::AbstractVector)
    current_best_learning(examples::AbstractVector, h::AbstractVector)

Apply the current-best-hypothesis learning algorithm (Fig. 19.2) on the given examples 'examples'
and hypothesis 'h' (an array of dictionaries where each Dict represents a disjunction). Return
a consistent hypothesis if possible, otherwise 'nothing' on failure.
"""
function current_best_learning(examples::AbstractVector, h::AbstractVector, prior_examples::AbstractVector)
    if (length(examples) == 0)
        return h;
    end

    local example::Dict = examples[1];

    push!(prior_examples, example);

    if (example_is_consistent(example, h))
        return current_best_learning(examples[2:end], h, prior_examples);
    elseif (example_is_false_positive(example, h))
        for h_prime in specializations(prior_examples, h)
            h_prime_prime = current_best_learning(examples[2:end], h_prime, prior_examples);
            if (!(typeof(h_prime_prime) <: Nothing))
                return h_prime_prime;
            end
        end
    elseif (example_is_false_negative(example, h))
        for h_prime in generalizations(prior_examples, h)
            h_prime_prime = current_best_learning(examples[2:end], h_prime, prior_examples);
            if (!(typeof(h_prime_prime) <: Nothing))
                return h_prime_prime;
            end
        end
    end
    return nothing;
end

function current_best_learning(examples::AbstractVector, h::AbstractVector)
    return current_best_learning(examples, h, []);
end

function version_space_update(V::AbstractVector, e::Dict)
    return collect(h for h in V if (example_is_consistent(e, h)));
end

function values_table(examples::AbstractVector)
    local values::Dict = Dict();
    for example in examples
        for (k, v) in example
            if (k == "GOAL")
                continue
            end
            local modifier::String = "!";
            if (example["GOAL"])
                modifier = "";
            end

            local modified_value::String = modifier * v;
            if (!(modified_value in get!(values, k, [])))
                push!(get!(values, k, []), modified_value);
            end
        end
    end
    return values;
end

function build_attribute_combinations(subset::Tuple, values::Dict)
    local h::AbstractVector = [];
    if (length(subset) == 1)
        k = values[subset[1]]
        h = collect([Dict([Pair(subset[1], v)])] for v in values[subset[1]]);
        return h;
    end
    for (i, attribute) in enumerate(subset)
        local rest::AbstractVector = build_attribute_combinations(subset[2:end], values);
        for value in values[attribute]
            local combination::Dict = Dict([Pair(attribute, value)]);
            for rest_item in rest
                local combination_prime::Dict = copy(combination);
                for dictionary in rest_item
                    merge!(combination_prime, dictionary);
                end
                push!(h, [combination_prime]);
            end
        end
    end
    return h;
end

function build_h_combinations(hypotheses::AbstractVector)
    local h::AbstractVector = [];
    local h_powerset::Set = setdiff!(generate_powerset(collect(1:length(hypotheses))), Set([()]));

    for subset in h_powerset
        local combination::AbstractVector = [];
        for index in subset
            append!(combination, hypotheses[index]);
        end
        push!(h, combination);
    end

    return h;
end

function all_hypotheses(examples::AbstractVector)
    local values::Dict = values_table(examples);
    local h_powerset::Set = setdiff!(generate_powerset(collect(keys(values))), Set([()]));
    local hypotheses::AbstractVector = [];
    for subset in h_powerset
        append!(hypotheses, build_attribute_combinations(subset, values));
    end
    append!(hypotheses, build_h_combinations(hypotheses));
    return hypotheses;
end

"""
    version_space_learning(examples::AbstractVector)

Return a version space for the given 'examples' by using the version space learning
algorithm (Fig. 19.3).
"""
function version_space_learning(examples::AbstractVector)
    local V::AbstractVector = all_hypotheses(examples);
    for example in examples
        if (length(V) != 0)
            V = version_space_update(V, example);
        end
    end
    return V;
end

function is_consistent_determination(A::AbstractVector, E::AbstractVector)
    local H::Dict = Dict();

    for example in E
        local attribute_values::Tuple = Tuple((collect(example[attribute] for attribute in A)...,));
        if (haskey(H, attribute_values))
            if (H[attribute_values] != example["GOAL"])
                return false;
            end
        end
        H[attribute_values] = example["GOAL"];
    end

    return true;
end

"""
    minimal_consistent_determination(E::AbstractVector, A::Set)

Return a set of attributes by using the algorithm for finding a minimal consistent
determination (Fig. 19.8).
"""
function minimal_consistent_determination(E::AbstractVector, A::Set)
    local n::Int64 = length(A);
    for i in 0:n
        for A_i in combinations(A, i);
            if (is_consistent_determination(A_i, E))
                return Set(A_i);
            end
        end
    end
    return nothing;
end

#=

    FOILKnowledgeBase is a knowledge base that consists of first order logic definite clauses,

    constant symbols, and predicate symbols used by foil().

=#
mutable struct FOILKnowledgeBase <: AbstractKnowledgeBase
    fol_kb::FirstOrderLogicKnowledgeBase
    constant_symbols::Set
    predicate_symbols::Set

    function FOILKnowledgeBase()
        return new(FirstOrderLogicKnowledgeBase(), Set(), Set());
    end

    function FOILKnowledgeBase(initial_clauses::Array{Expression, 1})
        local fkb::FOILKnowledgeBase = new(FirstOrderLogicKnowledgeBase(), Set(), Set());
        for clause in initial_clauses
            tell(fkb, clause);
        end
        return fkb;
    end
end

function tell(fkb::FOILKnowledgeBase, e::Expression)
    if (!is_logic_definite_clause(e))
        error("tell(): ", repr(e), " , is not a definite clause!");
    end

    tell(fkb.fol_kb, e);
    fkb.constant_symbols = union(fkb.constant_symbols, constant_symbols(e));
    fkb.predicate_symbols = union(fkb.predicate_symbols, predicate_symbols(e));

    return nothing;
end

function ask(fkb::FOILKnowledgeBase, e::Expression)
    return fol_bc_ask(fkb.fol_kb, e);
end

function retract(fkb::FOILKnowledgeBase, e::Expression)
    retract(fkb.fol_kb, e);
    nothing;
end

"""
    extend_example(fkb::FOILKnowledgeBase, example::Dict, literal::Expression)

Return an array of extended examples by extending the given example 'example' to satisfy
the given literal 'literal'.
"""
function extend_example(fkb::FOILKnowledgeBase, example::Dict, literal::Expression)
    local solution::AbstractVector = [];
    local substitutions::Tuple = ask(fkb, substitute(example, literal));
    for substitution in substitutions
        push!(solution, merge!(substitution, example));
    end
    return solution;
end

"""
    update_positive_examples(fkb::FOILKnowledgeBase, examples_positive::AbstractVector, extended_positive_examples::AbstractVector, target::Expression)

Return an array of uncovered positive examples given the positive examples 'positive_examples' and
the extended positive examples 'extended_positive_examples'.
"""
function update_positive_examples(fkb::FOILKnowledgeBase, examples_positive::AbstractVector, extended_positive_examples::AbstractVector, target::Expression)
    local uncovered_positive_examples::Array{Dict, 1} = Array{Dict, 1}();
    for example in examples_positive
        if (any((function(dict::Dict)
                    return all((dict[x] == example[x]) for x in keys(example));
                end),
                extended_positive_examples))
            tell(fkb, substitute(example, target));
        else
            push!(uncovered_positive_examples, example);
        end
    end
    return uncovered_positive_examples;
end

"""
    new_literals(fkb::FOILKnowledgeBase, clause::Tuple{Expression, AbstractVector})

Return a Tuple of literals given the known predicate symbols in the FOIL knowledge base 'fkb'
and the horn clause 'clause'.

Each literal in the returned literals share at least 1 variable with the given horn clause.
"""
function new_literals(fkb::FOILKnowledgeBase, clause::Tuple{Expression, AbstractVector})
    local share_known_variables::Set = variables(clause[1]);
    for literal in clause[2]
        union!(share_known_variables, variables(literal));
    end
    local result::Tuple = ();
    for (predicate, arity) in fkb.predicate_symbols
        local new_variables::Set = Set(collect(standardize_variables(expr("x"), standardize_variables_counter)
                                                for i in 1:(arity - 1)));
        for arguments in iterable_cartesian_product(fill(union(share_known_variables, new_variables), arity))
            if (any((variable in share_known_variables) for variable in arguments))
                result = Tuple((result..., Expression(predicate, arguments...,)));
            end
        end
    end
    return result;
end

"""
    choose_literal(fkb::FOILKnowledgeBase, literals::Tuple, examples::Tuple{AbstractVector, AbstractVector})

Return the best literal from the given literals 'literals' by comparing the information gained.
"""
function choose_literal(fkb::FOILKnowledgeBase, literals::Tuple, examples::Tuple{AbstractVector, AbstractVector})
    local information_gain::Function = (function(literal::Expression)
                                            local examples_positive::Int64 = length(examples[1]);
                                            local examples_negative::Int64 = length(examples[2]);
                                            local extended_examples::AbstractVector = collect(vcat(collect(extend_example(fkb, example, literal)
                                                                                for example in examples[i])...,)
                                                                                for i in 1:2);
                                            local extended_examples_positive::Int64 = length(extended_examples[1]);
                                            local extended_examples_negative::Int64 = length(extended_examples[2]);
                                            if ((examples_positive + examples_negative == 0) ||
                                                (extended_examples_positive + extended_examples_negative == 0))
                                                return (literal, -1);
                                            end
                                            local T::Int64 = 0;
                                            for example in examples[1]
                                                if (any((function(l_prime::Dict)
                                                            return all((l_prime[x] == example[x]) for x in keys(example));
                                                        end),
                                                        extended_examples[1]))
                                                    T = T + 1;
                                                end
                                            end
                                            return (literal, (T * log((extended_examples_positive * (examples_positive + examples_negative) + 0.0001)/((extended_examples_positive + extended_examples_negative) * examples_positive))));
                                        end);

    local gains::Tuple = map(information_gain, literals);
    return reduce((function(t1::Tuple, t2::Tuple)
                        if (getindex(t1, 2) < getindex(t2, 2))
                            return t2;
                        else
                            return t1;
                        end
                    end), gains)[1];
end

"""
    new_clause(fkb::FOILKnowledgeBase, examples::Tuple{AbstractVector, AbstractVector}, target::Expression)

Return a horn clause and the extended positive examples as Tuple.

The horn clause is represented as (consequent, array of antecendents).
"""
function new_clause(fkb::FOILKnowledgeBase, examples::Tuple{AbstractVector, AbstractVector}, target::Expression)
    local clause::Tuple = (target, Array{Expression, 1}());
    extended_examples = examples;
    while (length(extended_examples[2]) != 0)
        local literal::Expression = choose_literal(fkb, new_literals(fkb, clause), extended_examples);
        push!(clause[2], literal);
        extended_examples = (collect(vcat(collect(extend_example(fkb, example, literal)
                                                for example in extended_examples[i])...,)
                                    for i in 1:2)...,);
    end
    return (clause, extended_examples[1]);
end

"""
    foil(fkb::FOILKnowledgeBase, examples::Tuple{AbstractVector, AbstractVector}, target::Expression)

Return an array of horn clauses by using the FOIL algorithm (Fig. 19.12) on the given FOIL knowledge
base 'fkb', set of examples 'examples', and the target literal 'target'.
"""
function foil(fkb::FOILKnowledgeBase, examples::Tuple{AbstractVector, AbstractVector}, target::Expression)
    local clauses::AbstractVector = [];
    local positive_examples::AbstractVector;
    local negative_examples::AbstractVector;
    positive_examples, negative_examples = examples;

    while (length(positive_examples) != 0)
        local clause::Tuple;
        local positive_extended_examples::AbstractVector;
        clause, positive_extended_examples = new_clause(fkb, (positive_examples, negative_examples), target);
        # Remove postive examples covered by 'clause' from 'examples'
        positive_examples = update_positive_examples(fkb, positive_examples, positive_extended_examples, target);
        push!(clauses, clause);
    end
    return clauses;
end

max_stars_repo_name: mikhail-j/aima-julia

max_issues_repo_name: mikhail-j/aima-julia

max_forks_repo_name: mikhail-j/aima-julia


text:
"""
From MATLAB ODE Suite by Shampine
"""
@def rosenbrock2332unpack begin
  if typeof(cache) <: OrdinaryDiffEqMutableCache
    d = cache.tab.d
  else
    d = cache.d
  end
end

@def rosenbrock2332pre0 begin
  @rosenbrock2332unpack
  c1 = Θ*(1-Θ)/(1-2d)
  c2 = Θ*(Θ-2d)/(1-2d)
end

@muladd function _ode_interpolant(Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23ConstantCache,Rosenbrock23Cache,Rosenbrock32ConstantCache,Rosenbrock32Cache},idxs::Nothing,T::Type{Val{0}})
  @rosenbrock2332pre0
  @.. y₀ + dt*(c1*k[1] + c2*k[2])
end

@muladd function _ode_interpolant(Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23ConstantCache,Rosenbrock23Cache,Rosenbrock32ConstantCache,Rosenbrock32Cache},idxs,T::Type{Val{0}})
  @rosenbrock2332pre0
  @.. y₀[idxs] + dt*(c1*k[1][idxs] + c2*k[2][idxs])
end

@muladd function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23Cache,Rosenbrock32Cache},idxs::Nothing,T::Type{Val{0}})
  @rosenbrock2332pre0
  @.. out = y₀ + dt*(c1*k[1] + c2*k[2])
  out
end

@muladd function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23Cache,Rosenbrock32Cache},idxs,T::Type{Val{0}})
  @rosenbrock2332pre0
  @views @.. out = y₀[idxs] + dt*(c1*k[1][idxs] + c2*k[2][idxs])
  out
end

# First Derivative of the dense output
@def rosenbrock2332pre1 begin
  @rosenbrock2332unpack
  c1diff = (1-2*Θ)/(1-2*d)
  c2diff = (2*Θ-2*d)/(1-2*d)
end

@muladd function _ode_interpolant(Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23ConstantCache,Rosenbrock23Cache,Rosenbrock32ConstantCache,Rosenbrock32Cache},idxs::Nothing,T::Type{Val{1}})
  @rosenbrock2332pre1
  @.. c1diff*k[1] + c2diff*k[2]
end

@muladd function _ode_interpolant(Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23ConstantCache,Rosenbrock23Cache,Rosenbrock32ConstantCache,Rosenbrock32Cache},idxs,T::Type{Val{1}})
  @rosenbrock2332pre1
  @.. c1diff*k[1][idxs] + c2diff*k[2][idxs]
end

@muladd function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23Cache,Rosenbrock32Cache},idxs::Nothing,T::Type{Val{1}})
  @rosenbrock2332pre1
  @.. out = c1diff*k[1] + c2diff*k[2]
  out
end

@muladd function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Union{Rosenbrock23Cache,Rosenbrock32Cache},idxs,T::Type{Val{1}})
  @rosenbrock2332pre1
  @views @.. out = c1diff*k[1][idxs] + c2diff*k[2][idxs]
  out
end

"""
From MATLAB ODE Suite by Shampine
"""
@muladd function _ode_interpolant(Θ,dt,y₀,y₁,k,cache::Union{Rodas4ConstantCache,Rodas4Cache},idxs::Nothing,T::Type{Val{0}})
  Θ1 = 1 - Θ
  @.. Θ1*y₀ + Θ*(y₁ + Θ1*(k[1] + Θ*k[2]))
end

@muladd function _ode_interpolant(Θ,dt,y₀,y₁,k,cache::Union{Rodas4ConstantCache,Rodas4Cache},idxs,T::Type{Val{0}})
  Θ1 = 1 - Θ
  @.. Θ1*y₀[idxs] + Θ*(y₁[idxs] + Θ1*(k[1][idxs] + Θ*k[2][idxs]))
end

@muladd function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Rodas4Cache,idxs::Nothing,T::Type{Val{0}})
  Θ1 = 1 - Θ
  @.. out = Θ1*y₀ + Θ*(y₁ + Θ1*(k[1] + Θ*k[2]))
  out
end

@muladd function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Rodas4Cache,idxs,T::Type{Val{0}})
  Θ1 = 1 - Θ
  @views @.. out = Θ1*y₀[idxs] + Θ*(y₁[idxs] + Θ1*(k[1][idxs] + Θ*k[2][idxs]))
  out
end


text:
include("MackeyGlass.jl")

X, T = MGGenerator();

max_stars_repo_name: JuliaTagBot/MackeyGlassTimeSeries

max_issues_repo_name: JuliaTagBot/MackeyGlassTimeSeries

max_forks_repo_name: JuliaTagBot/MackeyGlassTimeSeries


text:
# K3S
#
# [-Q - ρI   Aᵀ      I           -I     ][ Δx  ]   [       -rc       ]
# [   A      δI      0            0     ][ Δy  ]   [       -rb       ]
# [   I      0   S_l⁻¹(X-L)       0     ][Δs_l2] = [-(x-l) + σμS_l⁻¹e]
# [  -I      0       0       S_u⁻¹(U-X) ][Δs_u2]   [-(u-x) + σμS_u⁻¹e]
export K3SStructuredParams

"""
Type to use the K3S formulation with a Krylov method, using the package
[`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
The outer constructor

    K3SStructuredParams(; uplo = :U, kmethod = :trimr, rhs_scale = true,
                         atol0 = 1.0e-4, rtol0 = 1.0e-4,
                         atol_min = 1.0e-10, rtol_min = 1.0e-10,
                         ρ0 =  sqrt(eps()) * 1e3, δ0 = sqrt(eps()) * 1e4,
                         ρ_min = 1e4 * sqrt(eps()), δ_min = 1e4 * sqrt(eps()),
                         mem = 20)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
The available methods are:
- `:tricg`
- `:trimr`
- `:gpmr`

The `mem` argument sould be used only with `gpmr`.
"""
mutable struct K3SStructuredParams <: NewtonParams
  uplo::Symbol
  kmethod::Symbol
  rhs_scale::Bool
  atol0::Float64
  rtol0::Float64
  atol_min::Float64
  rtol_min::Float64
  ρ0::Float64
  δ0::Float64
  ρ_min::Float64
  δ_min::Float64
  mem::Int
end

function K3SStructuredParams(;
  uplo::Symbol = :U,
  kmethod::Symbol = :trimr,
  rhs_scale::Bool = true,
  atol0::T = 1.0e-4,
  rtol0::T = 1.0e-4,
  atol_min::T = 1.0e-10,
  rtol_min::T = 1.0e-10,
  ρ0::T = sqrt(eps()) * 1e3,
  δ0::T = sqrt(eps()) * 1e4,
  ρ_min::T = 1e4 * sqrt(eps()),
  δ_min::T = 1e4 * sqrt(eps()),
  mem::Int = 20,
) where {T <: Real}
  return K3SStructuredParams(
    uplo,
    kmethod,
    rhs_scale,
    atol0,
    rtol0,
    atol_min,
    rtol_min,
    ρ0,
    δ0,
    ρ_min,
    δ_min,
    mem,
  )
end

mutable struct PreallocatedDataK3SStructured{
  T <: Real,
  S,
  L1 <: LinearOperator,
  L2 <: LinearOperator,
  L3 <: LinearOperator,
  Ksol <: KrylovSolver,
} <: PreallocatedDataNewtonStructured{T, S}
  rhs1::S
  rhs2::S
  rhs_scale::Bool
  regu::Regularization{T}
  δv::Vector{T}
  AI::L1
  Qreg::AbstractMatrix{T} # regularized Q
  QregF::LDLFactorizations.LDLFactorization{T, Int, Int, Int}
  Qregop::L2 # factorized matrix Qreg
  opBR::L3
  KS::Ksol
  atol::T
  rtol::T
  atol_min::T
  rtol_min::T
end

function opAIprod!(
  res::AbstractVector{T},
  nvar::Int,
  ncon::Int,
  ilow::Vector{Int},
  iupp::Vector{Int},
  nlow::Int,
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  A::AbstractMatrix{T},
  v::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  if β == 0
    res[(ncon + 1):(ncon + nlow)] .= @views α .* v[ilow]
    res[(ncon + nlow + 1):end] .= @views (-α) .* v[iupp]
  else
    res[(ncon + 1):(ncon + nlow)] .= @views α .* v[ilow] .+ β .* res[(ncon + 1):(ncon + nlow)]
    res[(ncon + nlow + 1):end] .= @views (-α) .* v[iupp] .+ β .* res[(ncon + nlow + 1):end]
  end
  if uplo == :U
    @views mul!(res[1:ncon], A', v, α, β)
  else
    @views mul!(res[1:ncon], A, v, α, β)
  end
end

function opAItprod!(
  res::AbstractVector{T},
  nvar::Int,
  ncon::Int,
  ilow::Vector{Int},
  iupp::Vector{Int},
  nlow::Int,
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  A::AbstractMatrix{T},
  v::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  if uplo == :U
    @views mul!(res, A, v[1:ncon], α, β)
  else
    @views mul!(res, A', v[1:ncon], α, β)
  end
  res[ilow] .+= @views α .* v[(ncon + 1):(ncon + nlow)]
  res[iupp] .-= @views α .* v[(ncon + nlow + 1):end]
end

function opBRK3Sprod!(
  res::AbstractVector{T},
  ncon::Int,
  nlow::Int,
  x_m_lvar::AbstractVector{T},
  uvar_m_x::AbstractVector{T},
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  δv::AbstractVector{T},
  v::AbstractVector{T},
  α::T,
  β::T,
) where {T <: Real}
  if β == zero(T)
    res[1:ncon] .= @views (α / δv[1]) .* v[1:ncon]
    res[(ncon + 1):(ncon + nlow)] .= @views α .* s_l ./ x_m_lvar .* v[(ncon + 1):(ncon + nlow)]
    res[(ncon + nlow + 1):end] .= @views α .* s_u ./ uvar_m_x .* v[(ncon + nlow + 1):end]
  else
    res[1:ncon] .= @views (α / δv[1]) .* v[1:ncon] .+ β .* res[1:ncon]
    res[(ncon + 1):(ncon + nlow)] .= @views α .* s_l ./ x_m_lvar .* v[(ncon + 1):(ncon + nlow)] .+
           β .* res[(ncon + 1):(ncon + nlow)]
    res[(ncon + nlow + 1):end] .=
      @views α .* s_u ./ uvar_m_x .* v[(ncon + nlow + 1):end] .+ β .* res[(ncon + nlow + 1):end]
  end
end

function update_kresiduals_historyK3S!(
  res::AbstractResiduals{T},
  Qreg::AbstractMatrix{T},
  AI::AbstractLinearOperator{T},
  δ::T,
  solx::AbstractVector{T},
  soly::AbstractVector{T},
  rhs1::AbstractVector{T},
  rhs2::AbstractVector{T},
  s_l::AbstractVector{T},
  s_u::AbstractVector{T},
  x_m_lvar::AbstractVector{T},
  uvar_m_x::AbstractVector{T},
  nvar::Int,
  ncon::Int,
  nlow::Int,
  ilow::Vector{Int},
  iupp::Vector{Int},
) where {T <: Real}
  if typeof(res) <: ResidualsHistory
    @views mul!(res.Kres[1:nvar], Symmetric(Qreg, :U), solx)
    @views mul!(res.Kres[1:nvar], AI', soly, one(T), one(T))
    @views mul!(res.Kres[(nvar + 1):end], AI, solx)
    res.Kres[(nvar + 1):(nvar + ncon)] .+= δ .* soly[1:ncon]
    res.Kres[(nvar + ncon + 1):(nvar + ncon + nlow)] .+=
      @views solx[ilow] .+ x_m_lvar .* soly[(ncon + 1):(ncon + nlow)] ./ s_l
    res.Kres[(nvar + ncon + nlow + 1):end] .+=
      @views .-solx[iupp] .+ uvar_m_x .* soly[(ncon + nlow + 1):end] ./ s_u
    res.Kres[1:nvar] .-= rhs1
    res.Kres[(nvar + 1):end] .-= rhs2
  end
end

function PreallocatedData(
  sp::K3SStructuredParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}

  # init Regularization values
  if iconf.mode == :mono
    regu = Regularization(T(sp.ρ0), T(sp.δ0), T(sp.ρ_min), T(sp.δ_min), :classic)
  else
    regu =
      Regularization(T(sp.ρ0), T(sp.δ0), T(sqrt(eps(T)) * 1e0), T(sqrt(eps(T)) * 1e0), :classic)
  end
  δv = [regu.δ] # put it in a Vector so that we can modify it without modifying opK2prod!
  # LinearOperator to compute [A  I  -I] v
  AI = LinearOperator(
    T,
    id.ncon + id.nlow + id.nupp,
    id.nvar,
    false,
    false,
    (res, v, α, β) -> opAIprod!(
      res,
      id.nvar,
      id.ncon,
      id.ilow,
      id.iupp,
      id.nlow,
      pt.s_l,
      pt.s_u,
      fd.A,
      v,
      α,
      β,
      fd.uplo,
    ),
    (res, v, α, β) -> opAItprod!(
      res,
      id.nvar,
      id.ncon,
      id.ilow,
      id.iupp,
      id.nlow,
      pt.s_l,
      pt.s_u,
      fd.A,
      v,
      α,
      β,
      fd.uplo,
    ),
  )

  Qreg = fd.Q + regu.ρ_min * I
  QregF = ldl(Symmetric(Qreg, :U))

  rhs1 = similar(fd.c, id.nvar)
  rhs2 = similar(fd.c, id.ncon + id.nlow + id.nupp)
  kstring = string(sp.kmethod)
  if sp.kmethod == :gpmr
    # operator to model the square root of the inverse of Q
    QregF.d .= sqrt.(QregF.d)
    Qregop = LinearOperator(
      T,
      id.nvar,
      id.nvar,
      false,
      false,
      (res, v) -> ld_div!(res, QregF, v),
      (res, v) -> dlt_div!(res, QregF, v),
    )
    # operator to model the square root of the inverse of the bottom right block of K3.5
    opBR = LinearOperator(
      T,
      id.ncon + id.nlow + id.nupp,
      id.ncon + id.nlow + id.nupp,
      true,
      true,
      (res, v, α, β) -> opsqrtBRK3Sprod!(
        res,
        id.ncon,
        id.nlow,
        itd.x_m_lvar,
        itd.uvar_m_x,
        pt.s_l,
        pt.s_u,
        δv,
        v,
        α,
        β,
      ),
    )
  else
    # operator to model the inverse of Q
    Qregop = LinearOperator(T, id.nvar, id.nvar, true, true, (res, v) -> ldiv!(res, QregF, v))
    # operator to model the inverse of the bottom right block of K3.5
    opBR = LinearOperator(
      T,
      id.ncon + id.nlow + id.nupp,
      id.ncon + id.nlow + id.nupp,
      true,
      true,
      (res, v, α, β) -> opBRK3Sprod!(
        res,
        id.ncon,
        id.nlow,
        itd.x_m_lvar,
        itd.uvar_m_x,
        pt.s_l,
        pt.s_u,
        δv,
        v,
        α,
        β,
      ),
    )
  end

  KS = init_Ksolver(AI', rhs1, sp)

  return PreallocatedDataK3SStructured(
    rhs1,
    rhs2,
    sp.rhs_scale,
    regu,
    δv,
    AI,
    Qreg,
    QregF,
    Qregop,
    opBR,
    KS,
    T(sp.atol0),
    T(sp.rtol0),
    T(sp.atol_min),
    T(sp.rtol_min),
  )
end

function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK3SStructured{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}
  Δs_l = itd.Δs_l
  Δs_u = itd.Δs_u
  pad.rhs1 .= @views step == :init ? fd.c : dd[1:(id.nvar)]
  if step == :init && all(pad.rhs1 .== zero(T))
    pad.rhs1 .= one(T)
  end
  pad.rhs2[1:(id.ncon)] .=
    @views (step == :init && all(dd[(id.nvar + 1):end] .== zero(T))) ? one(T) :
           dd[(id.nvar + 1):end]
  pad.rhs2[(id.ncon + 1):(id.ncon + id.nlow)] .= (step == :init) ? one(T) : Δs_l ./ pt.s_l
  pad.rhs2[(id.ncon + id.nlow + 1):end] .= (step == :init) ? one(T) : Δs_u ./ pt.s_u
  if pad.rhs_scale
    rhsNorm = sqrt(norm(pad.rhs1)^2 + norm(pad.rhs2)^2)
    pad.rhs1 ./= rhsNorm
    pad.rhs2 ./= rhsNorm
  end
  ksolve!(
    pad.KS,
    pad.AI',
    pad.rhs1,
    pad.rhs2,
    pad.Qregop,
    pad.opBR,
    verbose = 0,
    atol = pad.atol,
    rtol = pad.rtol,
  )
  update_kresiduals_historyK3S!(
    res,
    pad.Qreg,
    pad.AI,
    pad.regu.δ,
    pad.KS.x,
    pad.KS.y,
    pad.rhs1,
    pad.rhs2,
    pt.s_l,
    pt.s_u,
    itd.x_m_lvar,
    itd.uvar_m_x,
    id.nvar,
    id.ncon,
    id.nlow,
    id.ilow,
    id.iupp,
  )
  if pad.rhs_scale
    kunscale!(pad.KS.x, rhsNorm)
    kunscale!(pad.KS.y, rhsNorm)
  end

  dd[1:(id.nvar)] .= @views pad.KS.x
  dd[(id.nvar + 1):end] .= @views pad.KS.y[1:(id.ncon)]
  Δs_l .= @views pad.KS.y[(id.ncon + 1):(id.ncon + id.nlow)]
  Δs_u .= @views pad.KS.y[(id.ncon + id.nlow + 1):end]

  return 0
end

# function update_upper_Qreg!(Qreg, Q, ρ)
#   n = size(Qreg, 1)
#   nnzQ = nnz(Q)
#   if nnzQ > 0
#     for i=1:n
#       diag_idx = Qreg.colptr[i+1] - 1
#       ptrQ = Q.colptr[i+1] - 1
#       Qreg.nzval[diag_idx] = (ptrQ ≤ nnzQ && Q.rowval[ptrQ] == i) ? Q.nzval[ptrQ] + ρ : ρ
#     end
#   end
# end

function update_pad!(
  pad::PreallocatedDataK3SStructured{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end
  # if cnts.k == 4
  # update_upper_Qreg!(pad.Qreg, fd.Q, pad.regu.ρ)
  # pad.Qreg[diagind(pad.Qreg)] .= fd.Q[diagind(fd.Q)] .+ pad.regu.ρ
  # ldl_factorize!(Symmetric(pad.Qreg, :U), pad.QregF)
  # end

  if pad.atol > pad.atol_min
    pad.atol /= 10
  end
  if pad.rtol > pad.rtol_min
    pad.rtol /= 10
  end

  pad.δv[1] = pad.regu.δ

  return 0
end


text:
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    @warnpcfail precompile(Tuple{TaskThunk})
    @warnpcfail precompile(Tuple{typeof(wait_changed), String})
    @warnpcfail precompile(Tuple{typeof(watch_package), PkgId})
    @warnpcfail precompile(Tuple{typeof(watch_includes), Module, String})
    @warnpcfail precompile(Tuple{typeof(watch_manifest), String})
    @warnpcfail precompile(Tuple{typeof(revise_dir_queued), String})
    @warnpcfail precompile(Tuple{typeof(revise_file_queued), PkgData, String})
    @warnpcfail precompile(Tuple{typeof(init_watching), PkgData, Vector{String}})
    @warnpcfail precompile(Tuple{typeof(add_revise_deps)})
    @warnpcfail precompile(Tuple{typeof(watch_package_callback), PkgId})

    @warnpcfail precompile(Tuple{typeof(revise)})
    @warnpcfail precompile(Tuple{typeof(revise_first), Expr})
    @warnpcfail precompile(Tuple{typeof(includet), String})
    @warnpcfail precompile(Tuple{typeof(track), Module, String})
    # setindex! doesn't fully precompile, but it's still beneficial to do it
    # (it shaves off a bit of the time)
    # See https://github.com/JuliaLang/julia/pull/31466
    @warnpcfail precompile(Tuple{typeof(setindex!), ExprsSigs, Nothing, RelocatableExpr})
    @warnpcfail precompile(Tuple{typeof(setindex!), ExprsSigs, Vector{Any}, RelocatableExpr})
    @warnpcfail precompile(Tuple{typeof(setindex!), ModuleExprsSigs, ExprsSigs, Module})
    @warnpcfail precompile(Tuple{typeof(setindex!), Dict{PkgId,PkgData}, PkgData, PkgId})
    @warnpcfail precompile(Tuple{Type{WatchList}})
    @warnpcfail precompile(Tuple{typeof(setindex!), Dict{String,WatchList}, WatchList, String})

    MI = CodeTrackingMethodInfo
    @warnpcfail precompile(Tuple{typeof(minimal_evaluation!), MI, Core.CodeInfo, Symbol})
    @warnpcfail precompile(Tuple{typeof(minimal_evaluation!), Any, MI, Core.CodeInfo, Symbol})
    @warnpcfail precompile(Tuple{typeof(methods_by_execution!), Any, MI, DocExprs, Module, Expr})
    @warnpcfail precompile(Tuple{typeof(methods_by_execution!), Any, MI, DocExprs, JuliaInterpreter.Frame, Vector{Bool}})
    @warnpcfail precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                             NamedTuple{(:mode,),Tuple{Symbol}},
                             typeof(methods_by_execution!), Function, MI, DocExprs, Module, Expr})
    @warnpcfail precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                             NamedTuple{(:skip_include,),Tuple{Bool}},
                             typeof(methods_by_execution!), Function, MI, DocExprs, Module, Expr})
    @warnpcfail precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                             NamedTuple{(:mode, :skip_include),Tuple{Symbol,Bool}},
                             typeof(methods_by_execution!), Function, MI, DocExprs, Module, Expr})
    @warnpcfail precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                             NamedTuple{(:mode,),Tuple{Symbol}},
                             typeof(methods_by_execution!), Function, MI, DocExprs, Frame, Vector{Bool}})
    @warnpcfail precompile(Tuple{typeof(Core.kwfunc(methods_by_execution!)),
                             NamedTuple{(:mode, :skip_include),Tuple{Symbol,Bool}},
                             typeof(methods_by_execution!), Function, MI, DocExprs, Frame, Vector{Bool}})

    mex = which(methods_by_execution!, (Function, MI, DocExprs, Module, Expr))
    mbody = bodymethod(mex)
    # use `typeof(pairs(NamedTuple()))` here since it actually differs between Julia versions
    @warnpcfail precompile(Tuple{mbody.sig.parameters[1], Symbol, Bool, Bool, typeof(pairs(NamedTuple())), typeof(methods_by_execution!), Any, MI, DocExprs, Module, Expr})
    @warnpcfail precompile(Tuple{mbody.sig.parameters[1], Symbol, Bool, Bool, Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:skip_include,),Tuple{Bool}}}, typeof(methods_by_execution!), Any, MI, DocExprs, Module, Expr})
    mfr = which(methods_by_execution!, (Function, MI, DocExprs, Frame, Vector{Bool}))
    mbody = bodymethod(mfr)
    @warnpcfail precompile(Tuple{mbody.sig.parameters[1], Symbol, Bool, typeof(methods_by_execution!), Any, MI, DocExprs, Frame, Vector{Bool}})

    @warnpcfail precompile(Tuple{typeof(hastrackedexpr), Expr})
    @warnpcfail precompile(Tuple{typeof(get_def), Method})
    @warnpcfail precompile(Tuple{typeof(parse_pkg_files), PkgId})
    if isdefined(Revise, :filter_valid_cachefiles)
        @warnpcfail precompile(Tuple{typeof(filter_valid_cachefiles), String, Vector{String}})
    end
    @warnpcfail precompile(Tuple{typeof(pkg_fileinfo), PkgId})
    @warnpcfail precompile(Tuple{typeof(push!), WatchList, Pair{String,PkgId}})
    @warnpcfail precompile(Tuple{typeof(pushex!), ExprsSigs, Expr})
    @warnpcfail precompile(Tuple{Type{ModuleExprsSigs}, Module})
    @warnpcfail precompile(Tuple{Type{FileInfo}, Module, String})
    @warnpcfail precompile(Tuple{Type{PkgData}, PkgId})
    @warnpcfail precompile(Tuple{typeof(Base._deleteat!), Vector{Tuple{Module,String,Float64}}, Vector{Int}})
    @warnpcfail precompile(Tuple{typeof(add_require), String, Module, String, String, Expr})
    @warnpcfail precompile(Tuple{Core.kwftype(typeof(maybe_add_includes_to_pkgdata!)),NamedTuple{(:eval_now,), Tuple{Bool}},typeof(maybe_add_includes_to_pkgdata!),PkgData,String,Vector{Pair{Module, String}}})

    for TT in (Tuple{Module,Expr}, Tuple{DataType,MethodSummary})
        @warnpcfail precompile(Tuple{Core.kwftype(typeof(Base.CoreLogging.handle_message)),NamedTuple{(:time, :deltainfo), Tuple{Float64, TT}},typeof(Base.CoreLogging.handle_message),ReviseLogger,LogLevel,String,Module,String,Symbol,String,Int})
    end
    return nothing
end


text:
# TODO: add more target distributions and make them iteratable

# Dimension of testing distribution
const D = 5
# Tolerance ratio
const TRATIO = Int == Int64 ? 1 : 2
# Deterministic tolerance
const DETATOL = 1e-3 * D * TRATIO
# Random tolerance
const RNDATOL = 5e-2 * D * TRATIO

using Distributions: logpdf, MvNormal, InverseGamma, Normal
using DiffResults: GradientResult, JacobianResult, value, gradient, jacobian
using ForwardDiff: gradient!, jacobian!
using Bijectors: link, invlink, logpdf_with_trans

const dm = zeros(D)
const dσ = ones(D)
const MVN = MvNormal(dm, dσ)
ℓπ(θ) = logpdf(MVN, θ)

# Manual implementation of gradient so that it support matrix mode
function ∂ℓπ∂θ(θ::AbstractVecOrMat)
    d = MVN
    diff = θ .- dm
    dims = θ isa AbstractVector ? (:) : 1
    v = -(D * log(2π) + 2 * sum(log.(dσ)) .+ sum(abs2, diff ./ dσ; dims=dims)) / 2
    v = θ isa AbstractMatrix ? vec(v) : v
    g = -diff
    return (v, g)
end

function ∂ℓπ∂θ_ad(θ::AbstractVector)
    res = GradientResult(θ)
    gradient!(res, ℓπ, θ)
    return (value(res), gradient(res))
end

function ∂ℓπ∂θ_ad(θ::AbstractMatrix)
    v = similar(θ, size(θ, 2))
    g = similar(θ)
    for i in 1:size(θ, 2)
        res = GradientResult(θ[:,i])
        gradient!(res, ℓπ, θ[:,i])
        v[i] = value(res)
        g[:,i] = gradient(res)
    end
    return (v, g)
end

function ∂ℓπ∂θ_viajacob(θ::AbstractMatrix)
    jacob = similar(θ)
    res = JacobianResult(similar(θ, size(θ, 2)), jacob)
    jacobian!(res, ℓπ, θ)
    jacob_full = jacobian(res)
    d, n = size(jacob)
    for i in 1:n
        jacob[:,i] = jacob_full[i,1+(i-1)*d:i*d]
    end
    return (value(res), jacob)
end

# For the Turing model
# @model gdemo() = begin
#     s ~ InverseGamma(2, 3)
#     m ~ Normal(0, sqrt(s))
#     1.5 ~ Normal(m, sqrt(s))
#     2.0 ~ Normal(m, sqrt(s))
#     return s, m
# end

function invlink_gdemo(θ)
    s = invlink(InverseGamma(2, 3), θ[1])
    m = θ[2]
    return [s, m]
end

function ℓπ_gdemo(θ)
    s, m = invlink_gdemo(θ)
    logprior = logpdf_with_trans(InverseGamma(2, 3), s, true) + logpdf(Normal(0, sqrt(s)), m)
    loglikelihood = logpdf(Normal(m, sqrt(s)), 1.5) + logpdf(Normal(m, sqrt(s)), 2.0)
    return logprior + loglikelihood
end

function ∂ℓπ∂θ_gdemo(θ)
    res = GradientResult(θ)
    gradient!(res, ℓπ_gdemo, θ)
    return (value(res), gradient(res))
end


text:
# ------------------------------------------------------------------------------------------
# # Condicionales
#
# En Julia, la sintaxis
#
# ```julia
# if *condición 1*
#     *opción 1*
# elseif *condición 2*
#     *opción 2*
# else
#     *opción 3*
# end
# ```
#
# Nos permite eventualmente evaluar una de nuestras opciones.
#

# Por ejemplo, tal vez queremos implementar la prueba de FizzBuzz: Dado un número N, imprime
# "Fizz" si N es divisible entre 3, "Buzz" si N es divisible entre 5, y "FizzBuzz" si N es
# divisible entre ambos 3 y 5. En cualquier otro caso, imprimo el número mismo.
# ------------------------------------------------------------------------------------------

N =

if (N % 3 == 0) & (N % 5 == 0)
    println("FizzBuzz")
elseif N % 3 == 0
    println("Fizz")
elseif N % 5 == 0
    println("Buzz")
else
    println(N)
end

# ------------------------------------------------------------------------------------------
# Ahora digamos que queremos regresar el mayor número de ambos. Escoge tus propios x y y
# ------------------------------------------------------------------------------------------

x =
y =

if x > y
    x
else
    y
end

# ------------------------------------------------------------------------------------------
# Para el último bloque, podemos usar el operador ternario, con la sintaxis
#
# ```julia
# a ? b : c
# ```
#
# que equivale a
#
# ```julia
# if a
#     b
# else
#     c
# end
# ```
# ------------------------------------------------------------------------------------------

(x > y) ? x : y

# ------------------------------------------------------------------------------------------
# Un truco relacionado es la evaluación de corto-circuito
#
# ```julia
# a && b
# ```
# ------------------------------------------------------------------------------------------

(x > y) && println(x)

(x < y) && println(y)

# ------------------------------------------------------------------------------------------
# Cuando escribimos `a && b`, `b` se ejecuta sólo si `a` se evalúa a `true`.
#

# Si `a` se evalúa a `false`, la expresión `a && b` regresa `false`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
# 5.1 Reescribe FizzBuzz sin usar `elseif`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 5.2 Reescribe FizzBuzz usando el operador ternario.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
import ..UncertainDatasets:ConstrainedUncertainDataset

"""
	constrain(udata::UncertainDataset,
		s::SamplingConstraint) -> ConstrainedUncertainDataset

Return a uncertain dataset by applying the constraint `s` to each
uncertain value in `udata`.
"""
constrain(udata::UncertainDataset, constraint::SamplingConstraint) =
	ConstrainedUncertainDataset([constrain(uval, constraint) for uval in udata])

"""
	constrain(udata::UncertainDataset,
		constraints::Vector{T}) where {T<:SamplingConstraint} -> ConstrainedUncertainDataset

Return a uncertain dataset by applying a different sampling constraint to each uncertain
value in `udata`.
"""
function constrain(udata::UncertainDataset, constraints::Vector{T}) where {T<:SamplingConstraint}
	if length(udata) != length(constraints)
		error("Number of sampling constraints must match length of dataset.")
	end

	n_vals = length(udata)

	ConstrainedUncertainDataset([constrain(udata[i], constraints[i]) for i in 1:n_vals])
end

"""
	constrain(udata::ConstrainedUncertainDataset,
		constraint::SamplingConstraint) -> ConstrainedUncertainDataset

Return a uncertain dataset by applying the constraint `s` to each
uncertain value in `udata`.
"""
constrain(udata::ConstrainedUncertainDataset, constraint::SamplingConstraint) =
	ConstrainedUncertainDataset([constrain(uval, constraint) for uval in udata])

"""
	constrain(udata::ConstrainedUncertainDataset,
		constraints::Vector{T}) where {T<:SamplingConstraint} -> ConstrainedUncertainDataset

Return a uncertain dataset by applying a different sampling constraint to each uncertain
value in `udata`.
"""
function constrain(udata::ConstrainedUncertainDataset,
		constraints::Vector{T}) where {T<:SamplingConstraint}

	if length(udata) != length(constraints)
		error("Number of sampling constraints must match length of dataset.")
	end

	n_vals = length(udata)

	ConstrainedUncertainDataset([constrain(udata[i], constraints[i]) for i in 1:n_vals])
end

export constrain


text:
using Parquet

path = "c:/data/Performance_2003Q3.txt.parquet"
#Parquet.metadata(path)
@time col = Parquet.read_column(path, 5);

for i in 1:31
    println(i)
    @time Parquet.read_column(path, i);
end

@time read_parquet(path);

path = "c:/git/parquet-data-collection/dsd50p.parquet"
@time adf = read_parquet(path);

@time adf = read_parquet(path, multithreaded=false);

using JDF: type_compress!

@time adf = type_compress!(DataFrame(read_parquet(path, multithreaded=false), copycols=false));

using Random: randstring
tbl = (
    int32 = rand(Int32, 1000),
    int64 = rand(Int64, 1000),
    float32 = rand(Float32, 1000),
    float64 = rand(Float64, 1000),
    bool = rand(Bool, 1000),
    string = [randstring(8) for i in 1:1000],
    int32m = rand([missing, rand(Int32, 10)...], 1000),
    int64m = rand([missing, rand(Int64, 10)...], 1000),
    float32m = rand([missing, rand(Float32, 10)...], 1000),
    float64m = rand([missing, rand(Float64, 10)...], 1000),
    boolm = rand([missing, true, false], 1000),
    stringm = rand([missing, "abc", "def", "ghi"], 1000)
);

tmpfile = tempname()*".parquet"

@time write_parquet(tmpfile, tbl);

path = tmpfile
@time adf = read_parquet(path);

all([all(c1 .=== c2) for (c1, c2) in zip(tbl, adf)])

using BenchmarkTools
@benchmark adf = read_parquet(path)

col_num = 1
@time col1 = Parquet.read_column(path, col_num);
col1

meta = Parquet.metadata(path);
par = ParFile(path);

nrows(par)

colnames(par)
close(par)

#@time tbl = Parquet.read_column.(Ref(path), 1:length(colnames(par)));

col1

col1[19:20]

last(col1)

uncompressed_data_io = col1[1]

encoded_data_header = Parquet._read_varint(uncompressed_data_io, UInt32)

using Debugger

filemetadata = Parquet.metadata(path);
Debugger.@enter Parquet.read_column(path, filemetadata, col_num);

col1
correct = getproperty(tbl, keys(tbl)[col_num])
all(ismissing.(col1) .== ismissing.(correct))
all(skipmissing(col1) .== skipmissing(correct))

using Test
using Base.Threads: @spawn

checkcol(path, n; multithreaded=true) = begin
    res = Vector{Any}(undef, n)
    if multithreaded
        for col_num in 1:n
            res[col_num] = @spawn Parquet.read_column(path, col_num);
        end
        return fetch.(res)
    else
        for col_num in 1:n
            println(col_num)
            res[col_num] = Parquet.read_column(path, col_num);
        end
        return res
    end
end

@time checkcol(path, 31, multithreaded=true);
@time checkcol(path, 31, multithreaded=false);

@time checkcol(path, 12, multithreaded=false);

using Base.Threads: @spawn
read1(path, n) = begin
    result = Vector{Any}(undef, length(n))
    for i in n
        result[i] = @spawn Parquet.read_column(path, i)
    end
    fetch.(result)
end

@time a = read1(path, 1:5)

using DataFrames

@time ba=DataFrame(a, copycols=false)
@time ba=DataFrame(a)

b1

import Base: add_int
@edit Base.add_int(100, 1)

add_int


text:
layout(::Layout, body) =
    """
<html>
<head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    <title>Hello Livecoding</title>
    <link rel="stylesheet" href="/css/style.css"/>
</head>
<body>
<div>
$body
</div>
</body>
</html>"""


text:
mutable struct GaussianFeature
    _means::AbstractArray{Float64,1} # the list of means
    _var::Float64 # the variance of gaussian
    _n_kernels::Int64 # cache the number of kernels
    function GaussianFeature(means::AbstractArray{Float64,1}, var::Float64)
        new(means, var, size(means)[1])
    end
end

function _gauss(mean::Float64, var::Float64, x::AbstractArray{Float64,1})
    return exp.(-(0.5 / var) .* (x .- mean) .^ 2)
end

function transform(feature::GaussianFeature, x::AbstractArray{Float64,1})
    # this feature transforms R^1 to R^1
    # x is an array of size n_samples
    # returns Phi = [phi_1(x),,,, phi_N(x)] where phi_i(x) is a gaussian centered around means[i]
    n_samples = size(x)[1]
    n_kernels = feature._n_kernels
    Phi = zeros(n_samples, n_kernels)
    for i = 1:n_kernels
        mean = feature._means[i]
        var = feature._var
        Phi[:, i] = _gauss(mean, var, x)
    end

    return Phi
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
function DiffusionDefinition.diffusion_message(::Val{:Prokaryote})
    println("* * * * * * * * * * ")
    println("Prokaryote diffusion")
    println("* * * * * * * * * * ")
    println("Description")
    println("------------")
    println("Chemical Langevin equation for a simple system describing ")
    println("production of a protein that is repressing its own production.")
    println("It is represented by a 4-dimensional diffusion driven by an ")
    println("8-dimensional Wiener process:")
    println("\ndXₜ=S[θ .× h(Xₜ)]dt + S ⊙ √. ( θ .× h(X_t) ) dWₜ,\n")
    println("where the custom operation ⊙:ℝᵈˣᵏ→ℝᵈˣᵏ is defined via:")
    println("\n(M ⊙ μ)ᵢⱼ = Mᵢⱼμⱼ, i=1,…,d; j=1,…,k,\n")
    println("S is the stoichiometry matrix:")
    println("""S=[
    0   0   1   0   0   0  -1   0
    0   0   0   1  -2   2   0  -1
   -1   1   0   0   1  -1   0   0
   -1   1   0   0   0   0   0   0
]"""
)
    println("and the function h is given by:")
    println("\nh(x) = (x₃x₄, K-x₄, x₄, x₁, x₂(x₂-1)/2, x₃, x₁, x₂)ᵀ.\n")
    println("The underlying process X represents:")
    println("X = (RNA, P, P₂, DNA)")
    println("The model has eight parameters:")
    println("\tc1 c2 c3 c4 c5 c6 c7 c8 K")
    println("The first eight parameters represent the rate of chemical reactions")
    println("The last represent the total number of the DNA strands")
    println("\nConstructors")
    println("------------")
    println("You may initialize the process either by passing the parameters")
    println("in the order above, as in:")
    println("c₁, c₂, c₃, c₄, c₅, c₆, c₇, c₈, K = ...")
    println("P = Prokaryote(c₁, c₂, c₃, c₄, c₅, c₆, c₇, c₈, K)")
    println("or by keyword arguments")
    println("p = (c1=..., c2=..., c3=..., c4=..., c5=..., c6=..., c7=..., c8=..., K=...)")
    println("P = Prokaryote(;p...)")
    println()
end

@diffusion_process Prokaryote{T} begin
    :dimensions
    process --> 4
    wiener --> 8

    :parameters
    c --> (8, T)
    K --> Float64

    :additional
    statespace --> LowerBoundedStateSpace((1,2,3,4), (0.0,1.0,0.0,0.0))
end

function bound_satisfied(P::Prokaryote, x)
    DiffusionDefinition._bound_satisfied(P, x) && (P.K > x[4])
end

@inline _aux₁(x, P::Prokaryote) = P.c5*x[2]*(x[2]-1)
@inline _aux₂(x, P::Prokaryote) = P.c6*x[3]
@inline _aux₃(x, P::Prokaryote) = P.c2*(P.K - x[4])
@inline _aux₄(x, P::Prokaryote) = P.c1*x[3]*x[4]
@inline _aux₁₂(x, P::Prokaryote) = 2.0*_aux₂(x, P) - _aux₁(x, P)
@inline _aux₃₄(x, P::Prokaryote) = _aux₃(x, P) - _aux₄(x, P)

function DiffusionDefinition.b(t, x, P::Prokaryote)
    k1 = _aux₁₂(x, P)
    k2 = _aux₃₄(x, P)
    @SVector [
        P.c3*x[4] - P.c7*x[1],
        P.c4*x[1] + k1 - P.c8*x[2],
        k2 - 0.5 * k1,
        k2
    ]
end

function _σ_prokaryote(x, P)
    k₁ = sqrt(0.5*_aux₁(x, P))
    k₂ = sqrt(_aux₂(x, P))
    k₃ = sqrt(_aux₃(x, P))
    k₄ = sqrt(_aux₄(x, P))
    _O = zero(x[1]) # check if this breaks things for gradients

    @SMatrix [
        _O  _O  sqrt(P.c3*x[4]) _O              _O      _O      -sqrt(P.c7*x[1])    _O;
        _O  _O  _O              sqrt(P.c4*x[1]) -2.0*k₁ 2.0*k₂  _O                  -sqrt(P.c8*x[2]);
        -k₄ k₃  _O              _O              k₁      -k₂     _O                  _O;
        -k₄ k₃  _O              _O              _O      _O      _O                  _O
    ]
end

function _a_prokaryote(x, P)
    k₁ = _aux₁(x, P) + 2.0*_aux₂(x, P)
    k₂ = _aux₃(x, P) + _aux₄(x, P)
    s₁ = P.c3*x[4]+P.c7*x[1]
    s₂ = P.c4*x[1]+P.c8*x[2]
    _O = zero(x[1]) # check if this breaks things for gradients

    @SMatrix [
        s₁  _O      _O          _O;
        _O  s₂+2*k₁ -k₁         _O;
        _O  -k₁     0.5*k₁+k₂   k₂;
        _O  _O      k₂          k₂
    ]
end

DiffusionDefinition.σ(t, x, P::Prokaryote) = _σ_prokaryote(x, P)
DiffusionDefinition.a(t, x, P::Prokaryote) = _a_prokaryote(x, P)


text:
#=
nest_module:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-06
=#

# main module
module NestModule
    name = "Julia"

    # sub module1
    module Hello
        import ..NestModule

        function hello(name)
            println("Hello $(name)!!")
        end
    end

    # sub module2
    module Bye
        import ..NestModule

        function byebye(name)
            println("Byebye $(name)!!")
        end
    end

    function main()
        println("NestModule.name = ", name)
        Hello.hello(name)
        Bye.byebye(name)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .NestModule
    NestModule.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
@testset "264.ugly-number-ii.jl" begin
    @test nth_ugly_number(10) == 12
    @test nth_ugly_number(153) == 6144
end