text:
import Gen
import HypothesisTests
import Images
import Colors

using HypothesisTests
using Images

@everywhere using Gen;

#The locations of the objects
ObjectsPositions = [[54., 5, -35], [54., 5, -10], [54., 5, 15], [54., 5, 40],
                    [18., 5, -35], [18., 5, -10], [18., 5, 15], [18., 5, 40],
                    [-18., 5, -35], [-18., 5, -10], [-18., 5, 15], [-18., 5, 40],
                    [-54., 5, -35], [-54., 5, -10], [-54., 5, 15], [-54., 5, 40]]

function linear_planner_constant(startPos, destPos, num_pnts)
    dir = destPos - startPos
    speed = 5

    pnts = []

    pnt = deepcopy(startPos)
    for p in 1:num_pnts
        pnt = dir / norm(dir) * speed + pnt
        if speed > norm(dir)
            pnt = destPos
        end
        push!(pnts, pnt)
    end
    return pnts
end

function linear_planner_exponential(startPos, destPos, num_pnts)
    #assume make it 1/4 of the way (not constant speed!!)
    pnts = []
    colinearPnt = deepcopy(startPos)
    for p in 1:num_pnts
        for i=1:length(colinearPnt)
            colinearPnt[i] = (3 * destPos[i] + colinearPnt[i]) / 4
        end
        push!(pnts, colinearPnt)
    end
    return pnts
end

function parabolic_planner_constant(startPos, destPos, num_pnts, a = .01, phirot = 0)
    xaxis = destPos - startPos
    xaxis /= norm(xaxis)

    thetay = atan2(xaxis[2], xaxis[1]) + pi/2 #make orthogonal
    phiy = acos(xaxis[3]) + phirot #rotate if desired

    yaxis = [cos(thetay)*sin(phiy), sin(thetay)*sin(phiy), sin(phiy)]
    yaxis /= norm(yaxis)

    zaxis = cross(xaxis, yaxis)
    zaxis /= norm(zaxis)

    np = norm(startPos)
    nd = norm(destPos - startPos)

    #actual length is less
    l_inc = (np * nd + nd ^ 3 / 3 + (1 - a) / 2 * nd^2) / num_pnts
    l = l_inc

    pnts = []
    atEnd = false
    for i in 1:num_pnts
        #rs = roots(BigFloat.([-l, np, (1 - a) / 2, 1/3]))
        k = i * nd / num_pnts
        #set = false

        #for r in rs
        #    if isreal(r)
        #        if Float64(r) > 0
        #            k = Float64(r)
        #            set = true
        #            break
        #        end
        #    end
        #end

        #if (!set || atEnd)
        #    pnt = destPos
        #    atEnd = true
        #else
            pnt = startPos + k * xaxis + a * k * (nd - k) * zaxis
        #end

        push!(pnts, pnt)
        l += l_inc
    end
    return pnts
end

#Define the probabilistic program
@program move_hand_to_destination(planner::Function, path_length = 50) begin

    mv_covar = [1. 0 0 ; 0 1. 0; 0 0 1.]

    startPos = @tag(mvnormal([0, 50., -55.], mv_covar), "startPos")

    destinationObject = @tag(uniform_discrete(1, 16), "destinationObject")
    destinationPos = @tag(mvnormal(ObjectsPositions[destinationObject], mv_covar), "destinationPos")

    #move along path (plus some noise)
    pnts = planner(startPos, destinationPos, path_length)
    for i=1:path_length
        nextPos = @tag(mvnormal(pnts[i], mv_covar), "pathPos$i")
    end
end

#Sample likely goals with Inference and simulations from the simple planner
function infer_planner_samples(planner::Function, num_samples::Int, sir, path_length = 50, prior_trace = nothing)

    if sir > 0
        #dictionary of the traces and their scores
        traces = Vector{Trace}(num_samples)
        scores = Vector{Float64}(num_samples)

        #collect a sample
        for i=1:num_samples

            #create the trace
            t = Trace()

            #add priors if know
            if (prior_trace != nothing)
                for k in keys(prior_trace.constraints)
                    constrain!(t, k, prior_trace.constraints[k])
                end
                for k in keys(prior_trace.interventions)
                    intervene!(t, k, prior_trace.interventions[k])
                end
                #for i in keys(prior_trace.proposals)
                #    propose!(t, i, prior_trace.proposals[i])
                #end
            end

            #generate the sample
            scores[i] = @generate(t, move_hand_to_destination(planner, path_length))

            #score and store
            traces[i] = t
        end
        weights = exp.(scores - logsumexp(scores))
        weights = weights / sum(weights)
        return traces[rand(Distributions.Categorical(weights))]
    else
        #Nested inference MH

        currentTrace = Trace()

        #add priors if know
        if (prior_trace != nothing)
            for k in keys(prior_trace.constraints)
                constrain!(currentTrace, k, prior_trace.constraints[k])
            end
            for k in keys(prior_trace.interventions)
                intervene!(currentTrace, k, prior_trace.interventions[k])
            end
            #for i in keys(prior_trace.proposals)
            #    propose!(t, i, prior_trace.proposals[i])
            #end
        end

        #generate the sample
        old_score = @generate(currentTrace, move_hand_to_destination(planner, path_length))
        i = 1
        #collect a sample
        while i < num_samples

            #create the new sample trace
            t = deepcopy(currentTrace)

            #generate the sample
            new_score = @generate(t, move_hand_to_destination(planner, path_length))
            #if log(rand()) < new_score - old_score || i == 0
                old_score = new_score
                currentTrace = t
                i += 1
            #end
        end
        return currentTrace
    end
end

function manhattanDistance(objID1, objID2)
    return abs((objID1 - 1) % 4 - (objID2 - 1) % 4) + abs(Int(ceil(objID1 / 4) - ceil(objID2 / 4)))
end

function getMatrixLocation(objID)
    return (Int(ceil(objID / 4)), ((objID - 1) % 4) + 1)
end

planners = [linear_planner_constant, parabolic_planner_constant]
sampledStepExperiments = [.1, .2, .5, 1]
inferenceAlgorithms = [[1, true], [32, true], [128, true], #[512, true],
                        [1, false], [32, false]]#, [256, false], [1024, false]]
#For each object, load sample linear file, take varying amount of steps from it and
#attempt to infer the goal node via the framework

distributionsOverSamples = Dict()
for planner in planners
    distributionsOverSteps = Dict()
    for steps in sampledStepExperiments
        distributionsOverInference = Dict()
        for size in inferenceAlgorithms
            distributionsOverInference[size] = zeros(4, 4)
        end
        distributionsOverSteps[steps] = distributionsOverInference
    end
    distributionsOverSamples[planner] = distributionsOverSteps
end

entropyOverSamples = Dict()
for planner in planners
    entropyOverSteps = Dict()
    for steps in sampledStepExperiments
        entropyOverInference = Dict()
        for size in inferenceAlgorithms
            entropyOverInference[size] = zeros(4, 4)
        end
        entropyOverSteps[steps] = entropyOverInference
    end
    entropyOverSamples[planner] = entropyOverSteps
end

timeOverSamples = Dict()
for planner in planners
    timeOverSteps = Dict()
    for steps in sampledStepExperiments
        timeOverInference = Dict()
        for size in inferenceAlgorithms
            timeOverInference[size] = zeros(100,1)
        end
        timeOverSteps[steps] = timeOverInference
    end
    timeOverSamples[planner] = timeOverSteps
end

distribution = Dict()
for object in 1:16
    distribution[object] = 0
end

#Task configuration
print("Task configuration:\n\t")
for obj in 1:16
    print(obj, "\t")
    if (obj % 4 == 0)
        print("\n\t")
    end
end

linearData = false
customObj = false

for planner in planners
    println("TESTING ON PLANNER: $planner")

    for object in 1:16

        for i in ARGS
            try
                object = parse(Int, i)
                customObj = true
                println("Testing on user specified object:")
            catch
            end
            if i == "-l"
                linearData = true
                println("Testing on linear data set")
            end
        end
        println("Testing on object: $object")

        file = nothing
        if linearData
            file = open(joinpath(@__DIR__, "csv/linear/$(object + 16).csv"), "r")
        else
            file = open(joinpath(@__DIR__, "csv/$(object).csv"), "r")
        end
        humanPoints = []
        line = readline(file)
        while line != ""
            pnt = float(split(line, ",")[1:3])
            if linearData
                shift = [20, -20, 20]
                pnt += shift
            end
            push!(humanPoints, pnt)
            line = readline(file)
        end

        #number of points to skip each frame
        #skip = Int(floor(length(humanPoints) / 50))
        skip = 1

        for sampledSteps in sampledStepExperiments

            println("\tTesting with $sampledSteps % observed points")

            #set the observed values
            prior_trace = Trace()
            intervene!(prior_trace, "startPos", humanPoints[1])
            for i in 1:Int(floor(sampledSteps * length(humanPoints)))
                constrain!(prior_trace, "pathPos$i", humanPoints[skip * i])
            end

            #now collect the distribution of goals for a given sample size
            for inferenceAlgorithm in inferenceAlgorithms
                println("\t\t(planner = $planner, goal = $object, % observations = $sampledSteps, Inference samples = $inferenceAlgorithm):")

                #create the distribution from 50 Inference runs
                for infIt in 1:100
                    tic()
                    distribution[value(infer_planner_samples(planner, inferenceAlgorithm..., length(humanPoints),
                    prior_trace), "destinationObject")] += 1
                    timeOverSamples[planner][sampledSteps][inferenceAlgorithm][infIt] = toq()
                end

                #print cumulative
                #print("\n\t")
                #for obj in 1:16
                #    print(distribution[obj] / 100, "\t")
                #    if (obj % 4 == 0)
                #        print("\n\t")
                #    end
                #end
                #print("\n")

                expectedDist = 0.0
                for obj in 1:16
                    expectedDist += distribution[obj] / 100 * manhattanDistance(object, obj)
                end
                distributionsOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(object)...] = expectedDist

                entropy = 0.0
                for obj in 1:16
                    if distribution[obj] > 0
                        entropy -= distribution[obj] / 100 * log(distribution[obj] / 100)
                    end
                end
                entropyOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(object)...] = entropy

                #reset distribution for next run
                for obj in 1:16
                    distribution[obj] = 0
                end
            end
        end

        if customObj
            break
        end
    end
end

#print expected distances
for planner in planners
    for sampledSteps in sampledStepExperiments
        for inferenceAlgorithm in inferenceAlgorithms
            println("\tExpected Distance for inference with $planner, $inferenceAlgorithm samples, $sampledSteps % observed steps):")

            print("Exp. Dist.:\n\t")
            max = 0
            for obj in 1:16
                print(distributionsOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(obj)...], "\t")
                if (obj % 4 == 0)
                    print("\n\t")
                end
                if (distributionsOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(obj)...] > max)
                    max = distributionsOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(obj)...]
                end
            end
            save("images/eerr-$planner-$(Int(100 * sampledSteps))-$inferenceAlgorithm.png", Gray.(ones(4, 4) - distributionsOverSamples[planner][sampledSteps][inferenceAlgorithm] / max))

            max = 0
            print("Entropy:\n\t")
            for obj in 1:16
                print(entropyOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(obj)...], "\t")
                if (obj % 4 == 0)
                    print("\n\t")
                end
                if (entropyOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(obj)...] > max)
                    max = entropyOverSamples[planner][sampledSteps][inferenceAlgorithm][getMatrixLocation(obj)...]
                end
            end

            save("images/ent-$planner-$(Int(100 * sampledSteps))-$inferenceAlgorithm.png", Gray.(ones(4, 4) - entropyOverSamples[planner][sampledSteps][inferenceAlgorithm] / max))

            print("Average Time:$(sum(timeOverSamples[planner][sampledSteps][inferenceAlgorithm][:] / 100)) sec\n\t")
        end
    end
end
println("")

function printSigTable(mat, disp, savefilename)
    print(disp)
    for i in 1:size(mat,1)
        print("\n\t")
        for j in 1:size(mat,2)
            print(mat[i, j],"\t")
        end
    end
    println("\n")
    writecsv(joinpath("/home/liam/Documents/RSI/Probcomp/reaching_goal_inference/data", savefilename), mat)
end

println("Signficance for different planners (Wilcoxon signed rank test):")
planSigsErrs = ones(length(planners), length(planners))
planSigsEnts = ones(length(planners), length(planners))
planSigsTimes = ones(length(planners), length(planners))

for i in 1:length(planners)
    for j in 1:length(planners)
        if i == j
            continue
        end
        planner1 = planners[i]
        planner2 = planners[j]

        errs1 = distributionsOverSamples[planner1][.5][[32,1]]
        errs2 = distributionsOverSamples[planner2][.5][[32,1]]
        planSigsErrs[i, j] = pvalue(SignedRankTest(reshape(errs1, length(errs1)), reshape(errs2, length(errs2))))

        ents1 = entropyOverSamples[planner1][.5][[32,1]]
        ents2 = entropyOverSamples[planner2][.5][[32,1]]
        planSigsEnts[i, j] = pvalue(SignedRankTest(reshape(ents1, length(ents1)), reshape(ents2, length(ents2))))

        times1 = timeOverSamples[planner1][.5][[32,1]]
        times2 = timeOverSamples[planner2][.5][[32,1]]
        planSigsTimes[i, j] = pvalue(SignedRankTest(reshape(times1, length(times1)), reshape(times2, length(times2))))
    end
end

printSigTable(planSigsErrs, "Expected Error Significance (wrt planner):", "sig/plansigerr.csv")
printSigTable(planSigsEnts, "Entropy Significance (wrt planner):", "sig/plansigent.csv")
printSigTable(planSigsTimes, "Time Significance (wrt planner):", "sig/plansigt.csv")

println("Signficance for varying observation sizes (Wilcoxon signed rank test):")
obsSigsErrs = ones(length(sampledStepExperiments), length(sampledStepExperiments))
obsSigsEnts = ones(length(sampledStepExperiments), length(sampledStepExperiments))
obsSigsTimes = ones(length(sampledStepExperiments), length(sampledStepExperiments))

for i in 1:length(sampledStepExperiments)
    for j in 1:length(sampledStepExperiments)
        if i == j
            continue
        end

        sampledSteps1 = sampledStepExperiments[i]
        sampledSteps2 = sampledStepExperiments[j]

        errs1 = distributionsOverSamples[linear_planner_constant][sampledSteps1][[32,1]]
        errs2 = distributionsOverSamples[linear_planner_constant][sampledSteps2][[32,1]]
        obsSigsErrs[i, j] = pvalue(SignedRankTest(reshape(errs1, length(errs1)), reshape(errs2, length(errs2))))
        ents1 = entropyOverSamples[linear_planner_constant][sampledSteps1][[32,1]]
        ents2 = entropyOverSamples[linear_planner_constant][sampledSteps2][[32,1]]
        obsSigsEnts[i, j] = pvalue(SignedRankTest(reshape(ents1, length(ents1)), reshape(ents2, length(ents2))))

        times1 = timeOverSamples[linear_planner_constant][sampledSteps1][[32,1]]
        times2 = timeOverSamples[linear_planner_constant][sampledSteps2][[32,1]]
        obsSigsTimes[i, j] = pvalue(SignedRankTest(reshape(times1, length(times1)), reshape(times2, length(times2))))
    end
end

printSigTable(obsSigsErrs, "Expected Error Significance (wrt observation sample size):", "sig/obssigerr.csv")
printSigTable(obsSigsEnts, "Entropy Significance (wrt observation sample size):", "sig/obssigent.csv")
printSigTable(obsSigsTimes, "Time Significance (wrt observation sample size):", "sig/obssigt.csv")

println("Signficance for varying Inference sample sizes (Wilcoxon signed rank test):")
infSigsErrs = ones(length(inferenceAlgorithms), length(inferenceAlgorithms))
infSigsEnts = ones(length(inferenceAlgorithms), length(inferenceAlgorithms))
infSigsTimes = ones(length(inferenceAlgorithms), length(inferenceAlgorithms))

for i in 1:length(inferenceAlgorithms)
    for j in 1:length(inferenceAlgorithms)
        if i == j
            continue
        end

        inferenceAlgorithm1 = inferenceAlgorithms[i]
        inferenceAlgorithm2 = inferenceAlgorithms[j]

        errs1 = distributionsOverSamples[linear_planner_constant][.5][inferenceAlgorithm1]
        errs2 = distributionsOverSamples[linear_planner_constant][.5][inferenceAlgorithm2]
        infSigsErrs[i, j] = pvalue(SignedRankTest(reshape(errs1, length(errs1)), reshape(errs2, length(errs2))))

        ents1 = entropyOverSamples[linear_planner_constant][.5][inferenceAlgorithm1]
        ents2 = entropyOverSamples[linear_planner_constant][.5][inferenceAlgorithm2]
        infSigsEnts[i, j] = pvalue(SignedRankTest(reshape(ents1, length(ents1)), reshape(ents2, length(ents2))))

        times1 = timeOverSamples[linear_planner_constant][.5][inferenceAlgorithm1]
        times2 = timeOverSamples[linear_planner_constant][.5][inferenceAlgorithm2]
        infSigsTimes[i, j] = pvalue(SignedRankTest(reshape(times1, length(times1)), reshape(times2, length(times2))))
    end
end

printSigTable(infSigsErrs, "Expected Error Significance (wrt Inference sample size):", "sig/infsigerr.csv")
printSigTable(infSigsEnts, "Entropy Significance (wrt Inference sample size):", "sig/infsigent.csv")
printSigTable(infSigsTimes, "Time Significance (wrt Inference sample size):", "sig/infsigt.csv")


text:
function NumericalPropagatorBuilder(arg0::Orbit, arg1::ODEIntegratorBuilder, arg2::PositionAngle, arg3::jdouble)
    return NumericalPropagatorBuilder((Orbit, ODEIntegratorBuilder, PositionAngle, jdouble), arg0, arg1, arg2, arg3)
end

function NumericalPropagatorBuilder(arg0::Orbit, arg1::ODEIntegratorBuilder, arg2::PositionAngle, arg3::jdouble, arg4::AttitudeProvider)
    return NumericalPropagatorBuilder((Orbit, ODEIntegratorBuilder, PositionAngle, jdouble, AttitudeProvider), arg0, arg1, arg2, arg3, arg4)
end

function add_force_model(obj::NumericalPropagatorBuilder, arg0::ForceModel)
    return jcall(obj, "addForceModel", void, (ForceModel,), arg0)
end

function build_kalman_model(obj::NumericalPropagatorBuilder, arg0::List, arg1::List, arg2::ParameterDriversList)
    return jcall(obj, "buildKalmanModel", KalmanModel, (List, List, ParameterDriversList), arg0, arg1, arg2)
end

function build_ls_model(obj::NumericalPropagatorBuilder, arg0::Vector{IntegratedPropagatorBuilder}, arg1::List, arg2::ParameterDriversList, arg3::ModelObserver)
    return jcall(obj, "buildLSModel", BatchLSModel, (Vector{IntegratedPropagatorBuilder}, List, ParameterDriversList, ModelObserver), arg0, arg1, arg2, arg3)
end

function build_propagator(obj::NumericalPropagatorBuilder, arg0::Vector{jdouble})
    return jcall(obj, "buildPropagator", NumericalPropagator, (Vector{jdouble},), arg0)
end

function copy(obj::NumericalPropagatorBuilder)
    return jcall(obj, "copy", NumericalPropagatorBuilder, ())
end

function get_all_force_models(obj::NumericalPropagatorBuilder)
    return jcall(obj, "getAllForceModels", List, ())
end

function get_integrator_builder(obj::NumericalPropagatorBuilder)
    return jcall(obj, "getIntegratorBuilder", ODEIntegratorBuilder, ())
end

function get_mass(obj::NumericalPropagatorBuilder)
    return jcall(obj, "getMass", jdouble, ())
end

function set_mass(obj::NumericalPropagatorBuilder, arg0::jdouble)
    return jcall(obj, "setMass", void, (jdouble,), arg0)
end

ext: jl


text:
module Obj

import LsoBase

"""
    Objective(f, g [, h; f_b, g_b])

    Creates an objective function object with:
      f:     w  ->   f(w)
      g:     w  -> ∇ f(w)
      h:     w  -> ∇²f(w)
    f_b: (w, b) ->   f(w) (stochastic with indices b)
    g_b: (w, b) -> ∇ f(w) (stochastic with indices b)
    dim: Number of evaluations needed in order to compute
         full function value / gradient, e.g., the number
         of examples

    If f_b is not provided, f is used for f_b.
    If g_b is not provided, g is used for g_b.
    If h is not provided, an identity matrix is used.
"""
type Objective
      f::Function
      g::Function
      h::Function
    f_b::Function
    g_b::Function
    dim::Int
    function Objective(f::Function, g::Function,
                       h::Function = (w -> eye(length(w)));
                       f_b::Function = (w, b) -> f(w),
                       g_b::Function = (w, b) -> g(w),
                       dim::Int = -1)
        new(f, g, h, f_b, g_b, dim)
    end
end

"""
    f(obj, w [, b])

    Returns the function value of the objective function obj at point w.
    Optionally, use the index batch b to evaluate the function against.

    This helper function should ease the access to objects of type Obj.Objective.
"""
function f(obj::Objective, w::Array{Float64,1}, b::Array{Int,1}=Int[])
  if length(b) == 0
    return obj.f(w)
  else
    return obj.f_b(w, b)
  end
end

"""
    g(obj, w [, b])

    Returns the gradient of the objective function obj at point w.
    Optionally, use the index batch b to evaluate the function against.

    This helper function should ease the access to objects of type Obj.Objective.
"""
function g(obj::Objective, w::Array{Float64,1}, b::Array{Int,1}=Int[])
  if length(b) == 0
    return obj.g(w)
  else
    return obj.g_b(w, b)
  end
end

"""
    randbatch(obj, batchsize)

    Returns a random index batch of the given size for the objective function.
"""
function randbatch(obj::Objective, batchsize::Int=1)
  if batchsize == 1
    if obj.dim > 1
      return [rand(1:obj.dim)]
    else
      return Int[]
    end
  else
    if obj.dim < batchsize
      batchsize = obj.dim
    end
    return randperm(obj.dim)[1:batchsize]
  end
end

# include ./obj/*
LsoBase.includedir(dirname(@__FILE__)*"/obj")

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
"""
    BernoulliLikelihood(l::AbstractLink=LogisticLink())

Bernoulli likelihood is to be used if we assume that the
uncertainity associated with the data follows a Bernoulli distribution.
The link `l` needs to transform the input `f` to the domain [0, 1]

```math
    p(y|f) = Bernoulli(y | l(f))
```
On calling, this would return a Bernoulli distribution with `l(f)` probability of `true`.
"""
struct BernoulliLikelihood{Tl<:AbstractLink}
    invlink::Tl
end

BernoulliLikelihood() = BernoulliLikelihood(LogisticLink())

(l::BernoulliLikelihood)(f::Real) = Bernoulli(l.invlink(f))

(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(Bernoulli.(l.invlink.(fs)))


text:
using LinearAlgebra, CUDA, Flux, FastSolversForWeightedTV, TestImages, Test, Statistics
CUDA.allowscalar(false)

# Geometry
# flag_gpu = true
flag_gpu = false
T = Float64
n2d = (256,256)
n3d = (256,256,256)
u2d   = randn(T,          n2d); flag_gpu && (u2d   = u2d |> gpu)
u2d_c = randn(Complex{T}, n2d); flag_gpu && (u2d_c = u2d_c |> gpu)
u3d   = randn(T,          n3d); flag_gpu && (u3d   = u3d |> gpu)
u3d_c = randn(Complex{T}, n3d); flag_gpu && (u3d_c = u3d_c |> gpu)
rtol = T(1e-6)

# Operators
η2d = structural_mean(u2d)
P2d = structural_weight(u2d; η=η2d); flag_gpu && (P2d = P2d |> gpu)
η2d_c = structural_mean(u2d_c)
P2d_c = structural_weight(u2d_c; η=η2d_c); flag_gpu && (P2d_c = P2d_c |> gpu)
η3d = structural_mean(u3d)
P3d = structural_weight(u3d; η=η3d); flag_gpu && (P3d = P3d |> gpu)
η3d_c = structural_mean(u3d_c)
P3d_c = structural_weight(u3d_c; η=η3d_c); flag_gpu && (P3d_c = P3d_c |> gpu)

# Adjoint test (2d)
u = randn(T, n2d[1]-1, n2d[2]-1, 2); flag_gpu && (u = u |> gpu)
v = randn(T, n2d[1]-1, n2d[2]-1, 2); flag_gpu && (v = v |> gpu)
a = dot(P2d*u, v)
b = dot(u, adjoint(P2d)*v)
@test a ≈ b rtol = rtol

# Adjoint test (2d_c)
u = randn(Complex{T}, n2d[1]-1, n2d[2]-1, 2); flag_gpu && (u = u |> gpu)
v = randn(Complex{T}, n2d[1]-1, n2d[2]-1, 2); flag_gpu && (v = v |> gpu)
a = dot(P2d_c*u, v)
b = dot(u, adjoint(P2d_c)*v)
@test a ≈ b rtol = rtol

# Adjoint test (3d)
u = randn(T, n3d[1]-1, n3d[2]-1, n3d[3]-1, 3); flag_gpu && (u = u |> gpu)
v = randn(T, n3d[1]-1, n3d[2]-1, n3d[3]-1, 3); flag_gpu && (v = v |> gpu)
a = dot(P3d*u, v)
b = dot(u, adjoint(P3d)*v)
@test a ≈ b rtol = rtol

# Adjoint test (3d_c)
u = randn(Complex{T}, n3d[1]-1, n3d[2]-1, n3d[3]-1, 3); flag_gpu && (u = u |> gpu)
v = randn(Complex{T}, n3d[1]-1, n3d[2]-1, n3d[3]-1, 3); flag_gpu && (v = v |> gpu)
a = dot(P3d_c*u, v)
b = dot(u, adjoint(P3d_c)*v)
@test a ≈ b rtol = rtol

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
using Distributions
using Random
using SparseArrays

using CautiousEngine

# Get the absolute path of the experiment dir
EXPERIMENT_DIR = @__DIR__
experiment_type = "verification-example"
global_exp_dir = "$EXPERIMENT_DIR/$experiment_type/global"
specification_file = "foo"

#======================================================================================================
0. Define the Unknown Function, Process Noise, and the Full Dynamics Function
======================================================================================================#
# Full dynamics function with three control modes
A1 = [0.8 0.5; 0. 0.5]
A2 = [0.5 0.; -0.5 0.8]
g1 = (x) -> A1*x
g2 = (x) -> A2*x
g3 = (x) -> [x[1]+(sin(x[2])); 0.85*sqrt(abs(x[1]))*x[2]]
unknown_mode_list = [g1, g2, g3]
mode_names = ["mode1", "mode2", "mode3"]
known_part = nothing

# Define the zero-mean process noise
σ_proc = 0.01
measurement_noise_dist = Normal(0., σ_proc)
process_noise_dist = nothing                                # Process noise is only used in synthesis, not verification

#======================================================================================================
1.  Perform GP regression and construct IMDP for each unknown mode.
======================================================================================================#
exp_dir = global_exp_dir
res_dirs = Array{String,1}()                        # Array that holds the single-mode result directories
X = Dict("x1" => [-2., 2.], "x2" => [-2., 2.])      # Defines the compact safe set
grid_sizes = Dict("x1" => 0.25, "x2" => 0.25)       # Specifies the discretization coarseness in each dimension
safety_dims = ["x1", "x2"]                          # Dimensions that have safety bounds
dependency_dims = Dict("x1" => [1, 1],              # Defines dependencies on other states
                       "x2" => [1, 1])
number_of_datapoints = 100                          # Number of data points per mode to generate
run_exps_flag = true                                # Set to true to rerun the single-mode constructions
m_opt = -1                                          # Unused
bound_type = "rkhs-tight"                           # GP bound type
epsilon = -1.                                       # Distance bound parameter, set to -1 for automatic determination

for (i, mode) in enumerate(unknown_mode_list)
    system = mode_names[i]                          # Get the mode name
    Lf_bound = 1.0                                  # Lipschitz constant bound for the mode over X
    random_seed = 11 + i                            # Random seed for the random twister
    eta = σ_proc                                    # η parameter set in a naiive manner

    @info system, number_of_datapoints, σ_proc, bound_type, Lf_bound, random_seed

    # Setup experiment parameter structures (defined in CautiousEngine.jl)
    system_params = SystemParameters(system, mode, known_part, measurement_noise_dist, process_noise_dist, Lf_bound, dependency_dims, 2, 2)
    data_params = DataParameters(number_of_datapoints, m_opt, bound_type, σ_proc, epsilon, eta, safety_dims, -1., -1.)
    experiment_params = ExperimentParameters(exp_dir, experiment_type, specification_file, X, grid_sizes, "foo", random_seed, system_params, data_params)

    # If constructing the single-mode systems, run it!
    if run_exps_flag
        CautiousEngine.end_to_end_transition_bounds(experiment_params, reuse_gps_flag=false, reuse_regions_flag=false)
    end

    # Construct the directory where results will be
    res_dir = experiment_params.experiment_directory
    @info "Mode results can be found in $res_dir"
    push!(res_dirs, res_dir)
end

#======================================================================================================
2.  Construct the IMDP Labels and call the appropriate verification solver wrapper.
======================================================================================================#
default_label = "safe"                      # Default label of states
unsafe_label = "!safe"                      # The label of unsafe states
labels_dict = Dict("safe" => [X])           # Dictionary with key if label and value of corresponding extents
# Wrapper function
lbl = (state_extent; unsafe=false) -> general_label_fcn(state_extent, default_label, unsafe_label, labels_dict, unsafe=unsafe)
switched_system_dir = "$exp_dir/m$number_of_datapoints"
imdp = CautiousEngine.create_imdp_from_result_dirs(res_dirs, switched_system_dir, label_fn=lbl)

k_vals = [1, 2, 3, 4, 5, -1]                # Verification horizon (-1 means unbounded)
for k in k_vals
    verification_result_mat = CautiousEngine.Globally(imdp, "safe", k, "$switched_system_dir/switched-system/imdp-$k.txt")
    CautiousEngine.save_legacy_mats(verification_result_mat, "$switched_system_dir/switched-system", k)
end
CautiousEngine.plot_2d_verification_results("$switched_system_dir/switched-system")


text:
function _get_raster(::Type{EarthEnv}, ::Type{LandCover}, layer::Integer, full::Bool=false)
    1 ≤ layer ≤ 12 || throw(ArgumentError("The layer must be between 1 and 12"))

    filetype = full ? "complete" : "partial"
    path = joinpath(SimpleSDMLayers._layers_assets_path, _rasterpath(EarthEnv), _rasterpath(LandCover), filetype)
    isdir(path) || mkpath(path)

    root = "https://data.earthenv.org/consensus_landcover/"
    stem = full ? "with_DISCover/consensus_full_class_$(layer).tif" :
        "without_DISCover/Consensus_reduced_class_$(layer).tif"
    filename = "landcover_$(filetype)_$(layer).tif"

    if !isfile(joinpath(path, filename))
        Downloads.download(root * stem, joinpath(path, filename))
    end

    return joinpath(path, filename)
end

function _get_raster(::Type{EarthEnv}, ::Type{HabitatHeterogeneity}, layer::Integer, resolution::Int64=25)
    1 ≤ layer ≤ 14 || throw(ArgumentError("The layer must be between 1 and 14"))

    path = joinpath(SimpleSDMLayers._layers_assets_path, _rasterpath(EarthEnv), _rasterpath(HabitatHeterogeneity), string(resolution))
    isdir(path) || mkpath(path)

    # Layer names
    layernames = ["cv", "evenness", "range", "shannon", "simpson", "std", "Contrast", "Correlation", "Dissimilarity", "Entropy", "Homogeneity", "Maximum", "Uniformity", "Variance"]

    # Get the correct value for the numerical precision
    res = 16
    if resolution == 1
        if layer in [7, 9, 14]
            res = 32
        end
    end
    if resolution == 5
        if layer in [1, 7, 9, 14]
            res = 32
        end
    end
    if resolution == 25
        if layer in [7, 9, 10, 14]
            res = 32
        end
    end

    root = "https://data.earthenv.org/habitat_heterogeneity/"
    stem = "$(resolution)km/$(layernames[layer])_01_05_$(resolution)km_uint$(res).tif"
    filename = "$(layernames[layer])_$(resolution)km.tif"

    if !isfile(joinpath(path, filename))
        Downloads.download(root * stem, joinpath(path, filename))
    end

    return joinpath(path, filename)
end

function _get_raster(::Type{EarthEnv}, ::Type{Topography}, layer::Integer, resolution::Integer=50, source::String="GMTED", aggregation::String="mean")
    1 ≤ layer ≤ 16 || throw(ArgumentError("The layer must be between 1 and 16"))

    _src = Dict("GMTED" => "GMTED")
    _agr = Dict("mean" => "mn", "median" => "md", "minimum" => "mi", "maximum" => "ma", "std" => "sd ")
    _lay = ["elevation", "slope", "aspectcosine", "aspectsine", "eastness", "northness", "roughness", "tpi", "tri", "vrm", "dx" ,"dxx", "dy", "dyy", "pcurv", "tcurv"]
    _res = Dict(1 => "1KM", 10 => "10KM", 5 => "5KM", 50 => "50KM", 100 => "100KM")
    _sfx = layer == 1 ? _agr[aggregation] : "md"

    path = joinpath(SimpleSDMLayers._layers_assets_path, _rasterpath(EarthEnv), _rasterpath(Topography), _res[resolution])
    isdir(path) || mkpath(path)

    root = "https://data.earthenv.org/topography/"
    stem = "$(_lay[layer])_$(_res[resolution])$(_agr[aggregation])_$(_src[source])$(_sfx).tif"
    filename = stem

    if !isfile(joinpath(path, filename))
        Downloads.download(root * stem, joinpath(path, filename))
    end

    return joinpath(path, filename)
end


text:
__precompile__()

module Vega

    using JSON, ColorBrewer, KernelDensity, NoveltyColors, StatsBase, Parameters, Missings
    import Base: print, show

    export VegaAxis, VegaData, VegaMark, VegaPadding, VegaScale,
           VegaTransform, VegaVisualization, VegaFormat, VegaScene,
           VegaDataRef, VegaValueRef,
           VegaMarkPropertySet, VegaMarkProperties, VegaLegend, VegaStream, VegaMarkFrom, VegaPredicate

    export tojs, print

    export barplot, choropleth, lineplot, scatterplot, areaplot, heatmap, piechart,
           histogram, popchart, waterfall, wordcloud, groupedbar, rugplot, boxplot,
           bubblechart, streamplot, stemleaf, asterplot, jointplot, dotplot, ribbonplot,
           horizon, punchcard

    export xlab!, ylab!, xlim!, ylim!, title!, legend!, text!

    export default_scales!, default_axes!, add_data!
    export colorscheme!, coord_flip!, stroke!, hline!, vline!, hover!, jitter!, layer

    #Create base color library
    colorpalettes = merge(ColorBrewer.colorSchemes, NoveltyColors.ColorDict)

    #Import helper code
    include("primitives/definetypes.jl")
    include("primitives/tojs.jl")
    include("render.jl")

    # Intermediates
    include("intermediates/modifiers.jl")
    include("intermediates/defaults.jl")
    include("intermediates/add_data.jl")
    include("intermediates/xlab_ylab.jl")
    include("intermediates/xlim_ylim.jl")
    include("intermediates/title_text.jl")
    include("intermediates/hline_vline.jl")
    include("intermediates/layer.jl")

    # Higher-level API
    include("derived/choropleth.jl")
    include("derived/barplot.jl")
    include("derived/bubblechart.jl")
    include("derived/boxplot.jl")
    include("derived/groupedbar.jl")
    include("derived/lineplot.jl")
    include("derived/scatterplot.jl")
    include("derived/areaplot.jl")
    include("derived/heatmap.jl")
    include("derived/piechart.jl")
    include("derived/histogram.jl")
    include("derived/popchart.jl")
    include("derived/ribbonplot.jl")
    include("derived/rugplot.jl")
    include("derived/stemleaf.jl")
    include("derived/streamplot.jl")
    include("derived/waterfall.jl")
    include("derived/wordcloud.jl")
    include("derived/asterplot.jl")
    include("derived/jointplot.jl")
    include("derived/dotplot.jl")
    include("derived/horizon.jl")
    include("derived/punchcard.jl")

end


text:
#Write a function that checks whether an element occurs in a list.

list = [1,2,45,"A","B"]

function occ(y::Any,x::Vector{Any})
    y ∈ x
end

occ(1,list)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# creating load/risk over time model

using Pkg; Pkg.activate(".")
using PowerModels
using DataFrames

case = parse_file(joinpath(@__DIR__,"case5_risk_sys1.m"))
periods = 24

fname = case["name"]

# write load sin_wave
open(joinpath(@__DIR__, "$(fname)_load.csv"),"w") do io
    print(io, "period")
    for (id,load) in case["load"]; print(io,",load_$id"); end
    println(io)

    for i in 1:periods
        print(io,"$i")
        for (id,load) in case["load"]
            print(io, ",$(load["pd"]*0.5+0.5*sin(2*pi*i/periods))")
        end
        println(io)
    end
end

# write risk cos wave
comp_list = String[]
for comp_type in ["bus","branch","gen","load"]
    append!(comp_list, ["$(comp_type)_$(id)" for id in keys(case[comp_type])])
end

open(joinpath(@__DIR__, "$(fname)_risk.csv"),"w") do io
    print(io, "period")
    for comp in comp_list ; print(io,",$comp"); end
    println(io)

    for i in 1:periods
        print(io,"$i")
        for comp in comp_list
            print(io, ",$(rand()+0.5*cos(2*pi*i/periods))")
        end
        println(io)
    end
end

max_stars_repo_name: WISPO-POP/PowerModelsWildfireAnnex

max_issues_repo_name: WISPO-POP/PowerModelsWildfireAnnex

max_forks_repo_name: WISPO-POP/PowerModelsWildfireAnnex


text:
##INVESTIGATE COST OF SLICE SAMPLING WITHIN RE-ABC
using RareEventABC

##GENERATE OBSERVATIONS
srand(1);
d = 25;
yobs = randn(d) .* 3.0;
a = RareEventABC.IIDNormalDetails(yobs, d);

##INITIALISE SMC TUNING CHOICES
nparticles = 200;
parameters = [3.0];
schedule = doSMC(a, [3.0], 1.35, nparticles=200, quiet=true).schedule

##SPECIALISED VERSIONS OF SMC SUBROUTINES FOR THIS INVESTIGATION
##Slice sampling outputting the number of iterations required
function slicesample_(u::Array{Float64, 1}, dosim::Function, epsilon::Float64; w=1.0)
    ##Initialisation
    p = length(u)
    v = randn(p)
    a = -rand()*w
    b = w + a
    nits = 0
    ##Main loop
    while true
        z = a + rand() * (b-a)
        uprop = RareEventABC.reflect(u+z*v)
        sim, dist = dosim(uprop) ##Profiling suggests this line doesn't cause type-stability problems
        nits += 1
        if dist <= epsilon
            return uprop, sim, z, nits
        end
        if z < 0
            a = z
        else
            b = z
        end
    end
end

##Update particles storing mean number of slice sampling iterations required
function updateParticles_!(a::RareEventABC.IIDNormalDetails, p::RareEventABC.IIDNormalParticles, θ::Array{Float64, 1}, mean_ss_its::Vector; adaptive_w = true)
    function dosim(u::Array{Float64, 1})
        y = RareEventABC.u2y(u, θ[1])
        d = RareEventABC.getDistance(a, y)
        return (y, d)
    end
    z = zeros(p.nparticles)
    s = zeros(p.nparticles)
    for i in 1:p.nparticles
        if adaptive_w
            p.u[i], p.y[i], z[i], s[i] = slicesample_(p.u[i], dosim, p.ε, w=min(1.0,2.0*p.maxz))
        else
            p.u[i], p.y[i], z[i], s[i] = slicesample_(p.u[i], dosim, p.ε)
        end
    end
    p.maxz = maximum(abs(z))
    push!(mean_ss_its, mean(s))
    return nothing
end

##SMC - as usual but call functions above
function doSMC_(a::RareEventABC.ApplicationDetails, parameters::Array{Float64, 1}, schedule::Array{Float64, 1}, mean_ss_its::Vector; nparticles=200, llBound=-Inf, d=0, quiet=false, adaptive_w = true)
    ##Initialise SMC, including sampling from prior
    particles::RareEventABC.Particles = RareEventABC.Particles(a, nparticles, parameters)
    particles.schedule = schedule
    particles.llest = lgamma(1+d/2)-(d/2)*log(pi)-d*log(schedule[end]) ##log(1/V) where V is Lebesgue measure of a radius schedule[end] and dimension d
    particles.iteration = 0 ##Number of *completed* iterations
    particles.maxz = 1.0 ##Override what's currently in RareEventABC!

    ##Main loop
    nits = length(schedule)
    for ε in schedule
        particles.ε = ε
        dists::Array{Float64, 1} = RareEventABC.getDistances(a, particles)
        acc = (dists .<= ε)
        particles.llest += log(mean(acc))
        particles.iteration += 1
        ##quiet || report(a, particles, schedule[end], d, false)
        if particles.llest <= llBound
            particles.reachedBound = true
            break
        end
        if particles.iteration < nits
            ##These steps can be skipped on the last iteration
            RareEventABC.resample!(particles, acc)
            updateParticles_!(a, particles, parameters, mean_ss_its, adaptive_w=adaptive_w)
        end
    end
    particles
end

##A DETAILED COMPARISON
S_adaptive = Float64[];
S_nonadaptive = Float64[];
srand(2);
doSMC_(a, [3.0], schedule, S_adaptive, nparticles=200, quiet=true, adaptive_w = true).llest
doSMC_(a, [3.0], schedule, S_nonadaptive, nparticles=200, quiet=true, adaptive_w = false).llest

using DataFrames
using Plots; pgfplots()
using StatPlots
df = vcat(DataFrame(iteration=1:length(schedule)-1, mean_its=S_adaptive, w="adaptive"),
          DataFrame(iteration=1:length(schedule)-1, mean_its=S_nonadaptive, w="nonadaptive"));
default(guidefont=font(25), tickfont=font(25), legendfont=font(25), size=(500,500));
p1 = scatter(df, :iteration, :mean_its, group=:w, marker=:auto,
             xlabel="SMC iteration",
             ylabel="Mean slice sampling iterations",
             ylims=(1.0,11.0)); ##Hack to prevent legend getting in the way
savefig(p1, "slice_cost.pdf")

##REPEATED COMPARISONS TO SEE HOW TIME COST AND LIKELIHOOD ESTIMATES DIFFER
srand(0);
ll_adaptive = @time Float64[doSMC_(a, [3.0], schedule, S_adaptive, nparticles=200, quiet=true, adaptive_w = true).llest for i in 1:50];
ll_nonadaptive = @time Float64[doSMC_(a, [3.0], schedule, S_nonadaptive, nparticles=200, quiet=true, adaptive_w = false).llest for i in 1:50];

using HypothesisTests
UnequalVarianceTTest(exp(ll_adaptive), exp(ll_nonadaptive))

std(ll_adaptive)
std(ll_nonadaptive)


text:
struct PlotlyError
    msg::String
end

struct PlotlyCredentials
    username::String
    api_key::String
end

mutable struct PlotlyConfig
    plotly_domain::String
    plotly_api_domain::String
    plotly_streaming_domain::String
    plotly_proxy_authorization::Bool
    plotly_ssl_verification::Bool
    sharing::String
    world_readable::Bool
    auto_open::Bool
    fileopt::Symbol
end

const DEFAULT_CONFIG = PlotlyConfig(
    "https://plot.ly",
    "https://api.plot.ly/v2",
    "stream.plot.ly",
    false,
    true,
    "public",
    true,
    true,
    :create
)

function Base.merge(config::PlotlyConfig, other::AbstractDict)
    PlotlyConfig(
        [
            get(other, string(name), getfield(config, name))
            for name in fieldnames(PlotlyConfig)
        ]...
    )
end

Base.show(io::IO, config::PlotlyConfig) = dump(IOContext(io, :limit => true), config)

function Base.Dict(config::PlotlyConfig)
    Dict(k => getfield(config, k) for k in fieldnames(PlotlyConfig))
end

"""
    signin(username::String, api_key::String, endpoints=nothing)

Define session credentials/endpoint configuration, where endpoint is a Dict
"""
function signin(
        username::String, api_key::String,
        endpoints::Union{Nothing,AbstractDict}=nothing
    )
    global plotlycredentials = PlotlyCredentials(username, api_key)
    set_credentials_file(plotlycredentials)

    # if endpoints are specified both the base and api domains must be
    # specified
    if endpoints != nothing
        if !haskey(endpoints, "plotly_domain") || !haskey(endpoints, "plotly_api_domain")
            error("You must specify both the `plotly_domain` and `plotly_api_domain`")
        end
        global plotlyconfig = merge(DEFAULT_CONFIG, endpoints)
    end
end

"""
    get_credentials()

Return the session credentials if defined --> otherwise use .credentials specs
"""
function get_credentials()
    if !isdefined(Plotly, :plotlycredentials)

        creds = merge(get_credentials_file(), get_credentials_env())

        try
            username = creds["username"]
            api_key = creds["api_key"]

            global plotlycredentials = PlotlyCredentials(username, api_key)

        catch
            error("Please 'signin(username, api_key)' before proceeding. See
            http://plot.ly/API for help!")
        end
    end

    # will persist for the remainder of the session
    return plotlycredentials
end

"""
    get_config()

Return the session configuration if defined --> otherwise use .config specs
"""
function get_config()
    if !isdefined(Plotly, :plotlyconfig)
        config = get_config_file()
        global plotlyconfig = merge(DEFAULT_CONFIG, config)
    end

    # will persist for the remainder of the session
    return plotlyconfig
end

set_credentials_file(pc::PlotlyCredentials) = set_credentials_file(Dict("username" => pc.username, "api_key" => pc.api_key))

"""
    set_credentials_file(input_creds::AbstractDict)

Save Plotly endpoint configuration as JSON key-value pairs in
userhome/.plotly/.credentials. This includes username and api_key.
"""
function set_credentials_file(input_creds::AbstractDict)
    credentials_folder = joinpath(homedir(), ".plotly")
    credentials_file = joinpath(credentials_folder, ".credentials")

    # check to see if dir/file exists --> if not, create it
    !isdir(credentials_folder) && mkdir(credentials_folder)

    prev_creds = get_credentials_file()
    creds = merge(prev_creds, input_creds)

    # write the json strings to the cred file
    open(credentials_file, "w") do creds_file
        write(creds_file, JSON.json(creds))
    end
end

"""
    set_config_file(input_config::AbstractDict)

Save Plotly endpoint configuration as JSON key-value pairs in
userhome/.plotly/.config. This includes the plotly_domain, and
plotly_api_domain.
"""
function set_config_file(input_config::AbstractDict)
    config_folder = joinpath(homedir(), ".plotly")
    config_file = joinpath(config_folder, ".config")

    # check to see if dir/file exists --> if not create it
    !isdir(config_folder) && mkdir(config_folder)

    prev_config = get_config_file()
    config = merge(prev_config, input_config)

    # write the json strings to the config file
    open(config_file, "w") do config_file
        write(config_file, JSON.json(config))
    end
end

set_config_file(;kw...) = set_config_file(Dict(kw...))

"""
    set_config_file(config::PlotlyConfig)

Set the values in the configuration file to match the values in config
"""
set_config_file(config::PlotlyConfig) = set_config_file(Dict(config))

"""
    get_credentials_file()

Load user credentials informaiton as a dict
"""
function get_credentials_file()
    cred_file = joinpath(homedir(), ".plotly", ".credentials")
    isfile(cred_file) ? JSON.parsefile(cred_file) : Dict()
end

function get_credentials_env()
    out = Dict()
    keymap = Dict(
        "PLOTLY_USERNAME" => "username",
        "PLOTLY_APIKEY" => "api_key",
    )
    for k in ["PLOTLY_USERNAME", "PLOTLY_APIKEY"]
        if haskey(ENV, k)
            out[keymap[k]] = ENV[k]
        end
    end
    out
end

"""
    get_config_file()

Load endpoint configuration as a Dict
"""
function get_config_file()
    config_file = joinpath(homedir(), ".plotly", ".config")
    isfile(config_file) ? JSON.parsefile(config_file) : Dict()
end

max_stars_repo_name: plotly/Plotly-Julia


text:
#
# private.jl -
#
# Private module with stuff hidden to the end-user.  All methods, types and
# constants defined here are made available with the prefix `Priv.`.
#

# This method is meant to be extended.
function finalize end

# faster than convert(String,buf[1:end-1])
newstring(buf::Vector{<:Union{Int8,UInt8}}) =
    unsafe_string(pointer(buf))

newstring(buf::Vector{<:Union{Int8,UInt8}}, len::Integer) =
    unsafe_string(pointer(buf), len)


text:
import ..Common.SHAPE_FILL_ALPHA
import ..Common.SHAPE_COLOR

# Simulation
@recipe function f(sim::Simulation{1}, subplots=[1,2,3])
    x = map(a->a.vec[1],sim.x)
    n = sqrt.(sim.ε)

    layout --> (length(subplots),1)
    legend --> false

	for i ∈ 1:2
		if !(typeof(sim.boundary.bls[i])<:noBL)
			for j ∈ subplots
				@series begin
					subplot := j
					seriestype --> :path
					alpha --> 0
					if typeof(sim.boundary.bls[i])<:PML
						fillcolor --> :blue
					elseif typeof(sim.boundary.bls[i])<:cPML
						fillcolor --> :red
					end
					fillrange --> 0
					fillalpha --> SHAPE_FILL_ALPHA
					([sim.boundary.bls[i].start,
					  sim.boundary.bls[i].stop,
					  sim.boundary.bls[i].stop,
					  sim.boundary.bls[i].start,
					  sim.boundary.bls[i].start],

					  [-1000, -1000, 1000, 1000, -1000])
				end
			end
		end
	end
	if length(subplots) > 0
	    @series begin
			title --> L"{\rm Scattering\ Structure}"
			ylabel --> L"{\rm real}(n)"
	        subplot := subplots[1]
			ylims --> (min(1,minimum(real(n))), max(1,maximum(real(n)))) .+ (-.1,.1)
	        x, real(n)
	    end
	end
	if length(subplots) > 1
	    @series begin
			ylabel --> L"{\rm imag}(n)"
	        subplot := subplots[2]
			ylims --> (-1,1).*max(1,maximum(abs.(imag(n)))) .+ (-.1,.1)
	        x, imag(n)
	    end
	end
	if length(subplots) > 2
	    @series begin
			xlabel --> L"x"
			ylabel --> L"{\rm pump\ F}"
			subplot := subplots[3]
			ylims --> (-1,1).*max(1,maximum(abs.(sim.F))) .+ (-.1,.1)
	        x, sim.F
	    end
	end
end

# Scalar Field
@recipe function f(e::ScalarField{1}, by::Function)
    x = map(a->a.vec[1],e.pos)

	upper_lim = maximum(abs∘by,e.val)
	if by ∈ (abs,abs2)
		lower_lim = zero(upper_lim)
	else
		lower_lim = -upper_lim
	end

	m = size(e,2)
    layout --> (m,1)
    legend --> false
	ylims --> (lower_lim,upper_lim)
	for μ ∈ 1:m
	    @series begin
			subplot --> μ
			m>1 ? title --> latexstring("{\\rm Field}\\ $μ") : nothing
	        x, by.(e(μ).val)
	    end
	end
end

@recipe function f(sim::Simulation{1}, e::ScalarField{1})
	legend --> false
	m = size(e,2)
	layout --> (3,1+m)

	@series begin
		sim, [1,1+1+m,1+2(1+m)]
	end
    for μ ∈ 1:m
		@series begin
			subplot --> μ+1
			e(μ), real
		end
		@series begin
			subplot --> μ+1 + m+1
			e(μ), imag
		end
		@series begin
			subplot --> μ+1 + 2(m+1)
			e(μ), abs2
		end
	end
end

# Electric Field
@recipe function f(e::ElectricField{1}, by::Function)
    x = map(a->a.vec[1],e.pos)
    perm = sortperm(x)
	x = x[perm]

	upper_lim = maximum(abs∘by,e.val)
	if by ∈ (abs,abs2)
		lower_lim = zero(upper_lim)
	else
		lower_lim = -upper_lim
	end

	m = size(e,2)
    layout --> (3,m)
    legend --> false
	ylims --> (lower_lim,upper_lim)
	for μ ∈ 1:m
	    @series begin
			title --> latexstring("{\\rm Field}\\ $μ")
	        ylabel --> L"E_x"
	        subplot := 0m+μ
	        x, by.(e(μ).x[perm])
	    end
	    @series begin
	        ylabel --> L"E_y"
	        subplot := 1m+μ
	        x, by.(e(μ).y[perm])
	    end
	    @series begin
			xlabel --> L"x"
	        ylabel --> L"E_z"
	        subplot := 2m+μ
	        x, by.(e(μ).z[perm])
	    end
	end
end

@recipe function f(sim::Simulation{1},e::ElectricField{1}, by::Function)
    x = map(a->a.vec[1],e.pos)
    perm = sortperm(x)
    x = x[perm]

	legend --> false

	m = 1+size(e,2)
	layout --> (3,m)

	n = sqrt.(sim.ε[perm])
	@series begin
		title --> L"{\rm Scattering\ Structure}"
		ylabel --> L"{\rm real}(n)"
        subplot := 0m+1
        x, real(n)
    end
    @series begin
		ylabel --> L"{\rm imag}(n)"
        subplot := 1m+1
        x, imag(n)
    end
    @series begin
		xlabel --> L"x"
		ylabel --> L"{\rm pump\ F}"
		subplot := 2m+1
        x, sim.F[perm]
    end

	upper_lim = maximum(abs∘by,e.val)
	if by ∈ (abs,abs2)
		lower_lim = zero(upper_lim)
	else
		lower_lim = -upper_lim
	end
	ylims --> (lower_lim,upper_lim)

	for μ ∈ 1:m-1
	    @series begin
			title --> latexstring("{\\rm Field}\\ $μ")
	        ylabel --> L"E_x"
	        subplot := 0m+μ+1
	        x, by.(e(μ).x[perm])
	    end
	    @series begin
	        ylabel --> L"E_y"
	        subplot := 1m+μ+1
	        x, by.(e(μ).y[perm])
	    end
	    @series begin
			xlabel --> L"x"
	        ylabel --> L"E_z"
	        subplot := 2m+μ+1
	        x, by.(e(μ).z[perm])
	    end

	end
end


text:
#Functions to test for strong convergence of the SDE solvers.
#The number corresponds to the number of the formulas in Kloden & Platen 1992 - Chapter 4
import SDEs:integrate

"""
Approximating the integral of a function using the right hand side only.
"""
function rint(t::AbstractVector, x::AbstractVector)
    return sum(diff(t).*x[2:end])
end

"""
return the rms for the difference between 2 single-varible function (of t)
"""
function diff_rms(t::AbstractVector, f1vec::AbstractVector, f2vec::AbstractVector)
    fdiff = [norm(f1-f2) for (f1, f2) in zip(f1vec, f2vec)]
    rms = sqrt(rint(t, fdiff)/(t[end]-t[1]))
    return rms
end

function integrate(prob::Problem, x0::AbstractVector, t::AbstractVector, solver::Solver, seed::Integer; params...)
    return integrate(a!(prob), b!(prob), x0, t, solver, mgen(prob, x0), seed; params...)
end

function integrate(prob::Problem, x0::AbstractVector, t::AbstractVector, pce::Union{PCE, SDEs.PCEproto}, seed::Integer; params...)
    return integrate(a!(prob), b!(prob), x0, t, pce, mgen(prob, x0), seed; bbprime! = bbp!(prob), params...)
end

import DifferentialEquations
DiffEq = DifferentialEquations

struct DEQEM <: Solver end
function integrate(prob::Problem, x0::AbstractVector, t::AbstractVector, diffeqeuler::DEQEM, seed::Integer)
    nrp = zeros(length(x0), mgen(prob, x0))
    prob = DiffEq.SDEProblem(a!(prob), b!(prob), x0, (t[1], t[end]); noise_rate_prototype=nrp)
    sol = DiffEq.solve(prob, dt=t[2]-t[1], adaptive=false; save_noise=true, seed=seed)
    return sol.u[1:length(t)], sol.W.u[1:length(t)]
end

function test(prob::Problem, x0::AbstractVector, t::AbstractVector, solver::Solver, seed::Integer; params...)
    x, W = integrate(prob, x0, t, solver, seed; params...)
    x_analytic = solution(prob, x0, t, W)
    return x, x_analytic, W
end

function check(prob::Problem, x0::AbstractVector, t::AbstractVector, solver::Solver, seed::Integer; params...)
    x, x_analytic, _ = test(prob, x0, t, solver, seed; params...)
    return diff_rms(t, x, x_analytic)
end

function check(prob::Problem, x0::AbstractVector, t::AbstractVector, solver::Solver, seedlist::AbstractArray; params...)
    return mean([check(prob, x0, t, solver, i; params...) for i in seedlist])
end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
module PlanarConvexHulls

export
    ConvexHull,
    SConvexHull,
    DConvexHull,
    CW,
    CCW,
    vertices,
    num_vertices,
    area,
    centroid,
    is_ordered_and_strongly_convex,
    closest_point,
    jarvis_march!,
    hrep,
    hrep!

using LinearAlgebra
using StaticArrays
using StaticArrays: arithmetic_closure
using DocStringExtensions

include("order.jl")
include("util.jl")
include("exceptions.jl")
include("core_types.jl")
include("area.jl")
include("convexity_test.jl")
include("centroid.jl")
include("membership.jl")
include("closest_point.jl")
include("jarvis_march.jl")
include("hrep.jl")

end # module


text:
dimvals(x::NamedTuple,i) = x.axes[i]

_dname(::NamedTuple{N},i) where N = N[i]
_dname(::Tuple,i) = Symbol("Dim_",i)
dimname(x::NamedTuple, i) = _dname(x.axes,i)

#Special method to get all dim names
dimnames(x::NamedTuple) = _dnames(x.axes)
_dnames(::NamedTuple{N}) where N = N
_dnames(x) = ntuple(i->Symbol("Dim_",i), length(x))

getdata(x::NamedTuple) = x.values

function yaxcreate(::Type{<:NamedTuple},data,dnames,dvals,atts)
  axes = NamedTuple{dnames}(dvals)
  (axes = axes, values = data)
end


text:
#=
    We use two sets of infixable comparatives.

    The conventional symbols {==,!=,<,<=,>,>=} are defined using
       the Arb C library's implementation of eq,ne,lt,le,gt,ge.

    The predecessor/successor symbols {≃,≄,≺,≼,≻,≽} are defined from
       Hend Dawood's non-strict total ordering for interval values.
       (q.v. Hend's Master's thesis:
        Interval Mathematics Foundations, Algebraic Structures, and Applications)

=#

for (op,cop) in ((:(==), :(arb_eq)), (:(!=), :(arb_ne)),
                 (:(<=), :(arb_le)), (:(>=), :(arb_ge)),
                 (:(<), :(arb_lt)),  (:(>), :(arb_gt))  )
  @eval begin
    function ($op)(a::T, b::T) where {T <: ArbFloat}
        return Bool(ccall(@libarb($cop), Cint, (Ref{ArbFloat}, Ref{ArbFloat}), a, b) )
    end
    ($op)(a::ArbFloat{P}, b::ArbFloat{Q}) where {P,Q} = ($op)(promote(a,b)...)
    ($op)(a::T, b::R) where {T <: ArbFloat,R <: Real} = ($op)(promote(a,b)...)
    ($op)(a::R, b::T) where {T <: ArbFloat,R <: Real} = ($op)(promote(a,b)...)
  end
end

function (≃)(a::T, b::T) where {T <: ArbFloat}
    return Bool(ccall(@libarb(arb_eq), Cint, (Ref{ArbFloat}, Ref{ArbFloat}), a, b))
end
simeq(a::T, b::T) where {T <: ArbFloat} = (≃)(a,b)

function (≄)(a::T, b::T) where {T <: ArbFloat}
    return !Bool(ccall(@libarb(arb_eq), Cint, (Ref{ArbFloat}, Ref{ArbFloat}), a, b))
end
nsime(a::T, b::T) where {T <: ArbFloat} = (≄)(a,b)

function (⪰)(a::T, b::T) where {T <: ArbFloat}
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo < blo) || ((alo == blo) & (ahi <= bhi))
end
succeq(a::T, b::T) where {T <: ArbFloat} = (⪰)(a,b)

function (≻)(a::T, b::T) where {T <: ArbFloat} # (a ≼ b) & (a ≄ b)
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo < blo) || ((alo == blo) & (ahi < bhi))
end
succ(a::T, b::T) where {T <: ArbFloat} = (≻)(a,b)

function (⪯)(a::T, b::T) where {T <: ArbFloat}
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo > blo) || ((alo == blo) & (ahi >= bhi))
end
preceq(a::T, b::T) where {T <: ArbFloat} = (⪯)(a,b)

function (≺)(a::T, b::T) where {T <: ArbFloat}
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo > blo) || ((alo == blo) & (ahi > bhi))
end
prec(a::T, b::T) where {T <: ArbFloat} = (≺)(a,b)

# for sorted ordering
isequal(a::T, b::T) where {T <: ArbFloat} = !(a != b)
isless(a::T, b::T) where {T <: ArbFloat} = succ(a,b)
isequal(a::Nothing, b::T) where {T <: ArbFloat} = false
isequal(a::T, b::Nothing) where {T <: ArbFloat} = false
isless(a::Nothing, b::T) where {T <: ArbFloat} = true
isless(a::T, b::Nothing) where {T <: ArbFloat} = true

function max(x::T, y::T) where {T <: ArbFloat}
    (isnan(x) || isnan(y)) && return x
    if isinf(x)
       isposinf(x) && return x
       isneginf(x) && return y
    elseif isinf(y)
       isposinf(y) && return y
       isneginf(y) && return x
    end
    return (x + y + abs(x - y))/2
end

function min(x::T, y::T) where {T <: ArbFloat}
    (isnan(x) || isnan(y)) && return x
    if isinf(x)
       isposinf(x) && return y
       isneginf(x) && return x
    elseif isinf(y)
       isposinf(y) && return x
       isneginf(y) && return y
    end
    return (x + y - abs(x - y))/2
end
#=
min{T<:ArbFloat}(a::T, b::T) = smartvalue(a) < smartvalue(b) ? a : b
max{T<:ArbFloat}(a::T, b::T) = smartvalue(b) < smartvalue(a) ? a : b

function min2{T<:ArbFloat}(x::T, y::T)
    return
        if donotoverlap(x,y)
            x < y ? x : y
        else
            xlo, xhi = bounds(x)
            ylo, yhi = bounds(y)
            lo,hi = min(xlo, ylo), min(xhi, yhi)
            md = (hi+lo)/2
            rd = (hi-lo)/2
            midpoint_radius(md, rd)
        end
end

function max2{T<:ArbFloat}(x::T, y::T)
    return
        if donotoverlap(x,y)
            return x > y ? x : y
        else
            xlo, xhi = bounds(x)
            ylo, yhi = bounds(y)
            lo,hi = max(xlo, ylo), max(xhi, yhi)
            md = (hi+lo)/2
            rd = (hi-lo)/2
            return midpoint_radius(md, rd)
        end
end
=#
#=
function max{T<:ArbFloat}(x::T, y::T)
    return ((x>=y) | !(y<x)) ? x : y
end
=#

function minmax(x::T, y::T) where {T <: ArbFloat}
   return min(x,y), max(x,y) # ((x<=y) | !(y>x)) ? (x,y) : (y,x)
end

# experimental ≖ ≗
(eq)(x::T, y::T) where {T <: ArbFloat} = !(x != y)
(eq)(x::ArbFloat{P}, y::ArbFloat{Q}) where {P,Q} = (eq)(promote(x,y)...)
(eq)(x::T1, y::T2) where {T1 <: ArbFloat,T2 <: Real} = (eq)(promote(x,y)...)
(eq)(x::T2, y::T1) where {T1 <: ArbFloat,T2 <: Real} = (eq)(promote(x,y)...)
(≗)(x::T, y::T) where {T <: ArbFloat} = !(x != y)
(≗)(x::ArbFloat{P}, y::ArbFloat{Q}) where {P,Q} = (≗)(promote(x,y)...)
(≗)(x::T1, y::T2) where {T1 <: ArbFloat,T2 <: Real} = (≗)(promote(x,y)...)
(≗)(x::T2, y::T1) where {T1 <: ArbFloat,T2 <: Real} = (≗)(promote(x,y)...)
(neq)(x::T, y::T) where {T <: ArbFloat} = donotoverlap(x, y)
(neq)(x::ArbFloat{P}, y::ArbFloat{Q}) where {P,Q} = (donotoverlap)(promote(x,y)...)
(neq)(x::T1, y::T2) where {T1 <: ArbFloat,T2 <: Real} = (donotoverlap)(promote(x,y)...)
(neq)(x::T2, y::T1) where {T1 <: ArbFloat,T2 <: Real} = (donotoverlap)(promote(x,y)...)


text:
module ModelGraphSolvers

# using .BendersSolver

# using .ADMMSolver

export dual_decomposition_solve, DDOptimizer
#export benders_solve, ADMM_solve

include("DualDecompostion/DualDecompositionSolver.jl")
using .DualDecompSolver

end # module


text:
using CSV, DataFrames#, Distributed

function rerun_all(dir,filename_to_run)
    # @sync begin
    for (root, dirs, files) in walkdir(dir)
        for file in files
            if (file==filename_to_run)
                fpath = joinpath(root,file)
                try
                    #p = addprocs(1)
                    #@spawnat p[1] :(using CSV, DataFrames;
                    include(fpath); println("Success")
                    #)
                catch
                    @warn "Could not rerun file at "*fpath
                end
            end
        end
    end
    # end
    # rmprocs(procs()[2:end])
end

rerun_all(".","phd_ch_5_figs.jl")


text:
@testset "pdf of quadrature" begin

    gaussian(x, μ, σ) = exp(-0.5((x-μ)/σ)^2) / (σ*sqrt(2π))

    vacuum_state = VacuumState(rep=StateMatrix)
    v_μ = QSB.π̂ₓ_μ([0], vacuum_state)[1]
    v_σ = real(sqrt.(QSB.π̂ₓ²_μ([0], vacuum_state) .- v_μ.^2))[1]

    v(x) = gaussian(x, v_μ, v_σ)

    @test all(
        q_pdf(
            VacuumState(rep=StateMatrix, dim=100),
            LinRange(0, 2π, 2),
            LinRange(-10, 10, 100),
        )[i, :] ≈ v.(LinRange(-10, 10, 100))
        for i in 1:2
    )
end