text:
# Partial trace

using Convex
using SCS

import Base.sign
import Convex.monotonicity
import Convex.conic_form!
import Convex.curvature

export ptrace, partialtrace

"""
Type to be used internally by Convex.jl
"""

type PartialTraceAtom <: AbstractExpr
	head::Symbol
	id_hash::UInt64
	children::Tuple{AbstractExpr}
	size::Tuple{Integer, Integer}
	sys::Integer
	dims::Vector

	function PartialTraceAtom(x::AbstractExpr, sys::Integer, dims::Vector)
		if x.size[1] ≠ x.size[2]
			error("Only square matrices are supported")
		end

		if ! (1 ≤ sys ≤ length(dims))
			error("Invalid system, should between 1 and ", length(dims))
		end

		if x.size[1] ≠ prod(dims)
			error("Dimension of system doesnt correspond to dimension of subsystems")
		end

		children = (x, )
		newsize = (round(Int, x.size[2]/dims[sys]), round(Int, x.size[1]/dims[sys]))

		return new(:partialTrace,
		hash(children),
		children,
		newsize,
		sys,
		dims)
	end
end

function sign(x::PartialTraceAtom)
	 return NoSign()
	 return sign(x.children[0])
 end

function curvature(x::PartialTraceAtom)
  return ConstVexity()
end

function monotonicity(x::PartialTraceAtom)
	return (NoMonotonicity(),)
end

function conic_form!(x::PartialTraceAtom, unique_conic_forms::UniqueConicForms)
	if !has_conic_form(unique_conic_forms, x)
		sys = x.sys
		dims = x.dims
		function entry(ρ, j::Integer)
			bra = speye(1)
			ket = speye(1)
			i_sys = 1
			for dim in dims
				if i_sys == sys
					vO = sparsevec([j], [1], dim);
					bra = kron(bra, vO')
					ket = kron(ket, vO)
				else
					bra = kron(bra, speye(dim))
					ket = kron(ket, speye(dim))
				end
				i_sys += 1
			end
			return bra * ρ * ket
		end

		objective = conic_form!(sum([entry(x.children[1], j) for j in 1:dims[sys]]), unique_conic_forms)
		cache_conic_form!(unique_conic_forms, x, objective)
	end
	return get_conic_form(unique_conic_forms, x)
end

function partialtrace(x, systems::Vector, dims::Vector)
	dims = copy(dims)
	out = x
	for sys in systems
		out = partialtrace(out, sys, dims)
		dims[sys] = 1
	end
	return out
end

# General function, works in Convex.jl and with normal matrices

""" `rhoOut = partialtrace(rhoIn, sysNum, dimVec)`

Computes the partial trace of *rhoIn*, tracing out system number *sysNum*, where the dimensions of the systems are specified in the vector *dimVec*.

Example: rhoOut = partialtrace(epr, 2, [2 2])

Traces out the second qubit of a 2x2 state, here the EPR pair.
"""

function partialtrace(x, sys::Integer, dims::Vector)

	# Define a helper function returning matrix entries
	function entry(rho, j::Integer)
		bra = speye(1)
		ket = speye(1)
		i_sys = 1
		for dim in dims
			if i_sys == sys
				vO = sparsevec([j], [1], dim)
				bra = kron(bra, vO')
				ket = kron(ket, vO)
			else
				bra = kron(bra, speye(dim))
				ket = kron(ket, speye(dim))
			end
			i_sys += 1
		end
		return bra * rho * ket
	end

	# Compute the partial trace by picking the appropriate entries
	return sum([entry(x, j) for j in 1:dims[sys]])
end

ptrace(x::AbstractExpr, sys::Integer, dim::Vector) = partialtrace(x, sys, dim)
ptrace(x::AbstractArray, sys::Integer, dim::Vector) = partialtrace(x, sys, dim)


text:
"""
    JobDefinition

Stores the job definition arn including the revision.
"""
@auto_hash_equals struct JobDefinition
    arn::AbstractString

    function JobDefinition(name::AbstractString; aws_config::AbstractAWSConfig=global_aws_config())
        if startswith(name, "arn:")
            new(name)
        else
            arn = job_definition_arn(name; aws_config=aws_config)
            arn === nothing && error("No job definition ARN found for $name")
            new(arn)
        end
    end
end

"""
    job_definition_arn(
        definition_name::AbstractString;
        image::AbstractString="",
        role::AbstractString="",
        aws_config::AbstractAWSConfig=global_aws_config(),
    ) -> Union{AbstractString, Nothing}

Looks up the ARN (Amazon Resource Name) for the latest job definition that can be reused.
Returns a JobDefinition with the ARN that can be reused or `nothing`.

A job definition can only be reused if:

1. status = ACTIVE
2. type = container
3. image = the current job's image
4. jobRoleArn = the current job's role
"""
function job_definition_arn(
    definition_name::AbstractString;
    image::AbstractString="",
    role::AbstractString="",
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    response = describe_job_definition(definition_name; aws_config=aws_config)
    if !isempty(response["jobDefinitions"])

        latest = first(response["jobDefinitions"])
        for definition in response["jobDefinitions"]
            if definition["status"] == "ACTIVE" && definition["revision"] > latest["revision"]
                latest = definition
            end
        end

        if (
            latest["status"] == "ACTIVE" &&
            latest["type"] == "container" &&
            (latest["containerProperties"]["image"] == image || isempty(image)) &&
            (latest["containerProperties"]["jobRoleArn"] == role || isempty(role))
        )
            info(
                logger,
                string(
                    "Found previously registered job definition: ",
                    "\"$(latest["jobDefinitionArn"])\"",
                )
            )
            return latest["jobDefinitionArn"]
        end
    end

    notice(
        logger,
        string(
            "Did not find a previously registered ACTIVE job definition for ",
            "\"$definition_name\".",
        )
    )
    return nothing
end

"""
    register(
        definition_name::AbstractString;
        role::AbstractString="",
        image::AbstractString="",
        vcpus::Integer=1,
        memory::Integer=1024,
        cmd::Cmd=``,
        region::AbstractString="",
        parameters::Dict{String,String}=Dict{String, String}(),
    ) -> JobDefinition

Registers a new job definition.
"""
function register(
    definition_name::AbstractString;
    image::AbstractString="",
    role::AbstractString="",
    type::AbstractString="container",
    vcpus::Integer=1,
    memory::Integer=1024,
    cmd::Cmd=``,
    parameters::Dict{String, String}=Dict{String, String}(),
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    debug(logger, "Registering job definition \"$definition_name\"")
    input = OrderedDict(
        "parameters" => parameters,
        "containerProperties" => OrderedDict(
            "image" => image,
            "vcpus" => vcpus,
            "memory" => memory,
            "command" => cmd.exec,
            "jobRoleArn" => role,
        ),
    )

    response = @mock Batch.register_job_definition(definition_name, type, input; aws_config=aws_config)
    definition = JobDefinition(response["jobDefinitionArn"]; aws_config=aws_config)
    info(logger, "Registered job definition \"$(definition.arn)\"")
    return definition
end

"""
    deregister(job::JobDefinition)

Deregisters an AWS Batch job.
"""
function deregister(definition::JobDefinition; aws_config::AbstractAWSConfig=global_aws_config())
    debug(logger, "Deregistering job definition \"$(definition.arn)\"")
    resp = @mock Batch.deregister_job_definition(definition.arn; aws_config=aws_config)
    info(logger, "Deregistered job definition \"$(definition.arn)\"")
end

"""
    isregistered(definition::JobDefinition; aws_config=global_aws_config()) -> Bool

Checks if a JobDefinition is registered.
"""
function isregistered(definition::JobDefinition; aws_config::AbstractAWSConfig=global_aws_config())
    j = describe(definition; aws_config=aws_config)
    return any(d -> d["status"] == "ACTIVE", get(j, "jobDefinitions", []))
end

"""
    list_job_definitions(;aws_config=global_aws_config())

Get a list of `JobDefinition` objects via `Batch.decsribe_job_definitions()`.
"""
function list_job_definitions(;aws_config::AbstractAWSConfig=global_aws_config())
    job_definitions = Batch.describe_job_definitions(; aws_config=aws_config)["jobDefinitions"]

    return [JobDefinition(jd["jobDefinitionArn"]) for jd in job_definitions]
end

"""
    describe(definition::JobDefinition; aws_config=global_aws_config()) -> Dict

Describes a job definition as a dictionary. Requires the IAM permissions
"batch:DescribeJobDefinitions".
"""
function describe(definition::JobDefinition; aws_config::AbstractAWSConfig=global_aws_config())
    describe_job_definition(definition; aws_config=aws_config)
end

function describe_job_definition(definition::JobDefinition;
                                 aws_config::AbstractAWSConfig=global_aws_config())
    describe_job_definition(definition.arn; aws_config=aws_config)
end
function describe_job_definition(definition::AbstractString;
                                 aws_config::AbstractAWSConfig=global_aws_config())
    query = if startswith(definition, "arn:")
        Dict("jobDefinitions" => [definition])
    else
        Dict("jobDefinitionName" => definition)
    end
    return @mock Batch.describe_job_definitions(query; aws_config=aws_config)
end

ext: jl


text:
"""
`count_edges(channel; edge="rising",initial_count=0,direction="up") -> task`

create a NIDAQ counter input channel
"""
function count_edges(channel::String;
        edge::AbstractString="rising", initial_count::Integer=0, direction::AbstractString="up")
    t = CITask()
    if edge ∉ ("rising", "falling")
        error("edge must either be \"rising\" or \"falling\"")
    end
    if direction ∉ ("up", "down")
        error("direction must either be \"up\" or \"down\"")
    end
    catch_error( CreateCICountEdgesChan(t.th,
        Ref(codeunits(channel),1),
        Ref(codeunits(""),1),
        edge == "rising" ? Val_Rising : Val_Falling,
        UInt32(initial_count),
        direction == "up" ? Val_CountUp : Val_CountDown) )
    t
end

#= broken
function measure_duty_cycle(channel::String;  units::AbstractString="seconds")
    t = CITask()
    if units == "seconds"
        ret = CreateCIPulseChanTime(t.th,
                Ref(codeunits(channel),1), Ref(codeunits(""),1),
                2.0, 1000.0,
                Val_Seconds)
    elseif units == "ticks"
        ret = CreateCIPulseChanTicks(t.th,
                Ref(codeunits(channel),1),
                Ref(codeunits(""),1),
                Ref(codeunits(""),1),
                2.0, 1000.0)
    else
        error("units must either be \"seconds\" or \"ticks\"")
    end
    catch_error(ret)
    t
end
=#

"""
`quadrature_input(channel; z_enable=true) -> task`

create a NIDAQ counter input channel
"""
function quadrature_input(channel::String; z_enable::Bool=true)
    t = CITask()
    ret = CreateCIAngEncoderChan(t.th,
            Ref(codeunits(channel),1),
            Ref(codeunits(""),1),
            Val_X4,
            reinterpret(Bool32, UInt32(z_enable)),
            0.0,
            Val_AHighBHigh,
            Val_Ticks,
            UInt32(1), 0.0,
            Ref(codeunits(""),1))
    ret>0 && @warn(error(ret))
    ret<0 && error(error(ret))
    t
end

"""
`line_to_line(channel; units="seconds",edge1="rising",edge2="rising") -> task`

create a NIDAQ counter input channel
"""
function line_to_line(channel::String;
        units::AbstractString="seconds", edge1::AbstractString="rising", edge2::AbstractString="rising")
    t = CITask()
    if units ∉ ("seconds", "ticks")
        error("units must either be \"seconds\" or \"ticks\"")
    end
    if edge1 ∉ ("rising", "falling")
        error("edge1 must either be \"rising\" or \"falling\"")
    end
    if edge2 ∉ ("rising", "falling")
        error("edge2 must either be \"rising\" or \"falling\"")
    end
    catch_error( CreateCITwoEdgeSepChan(t.th,
            Ref(codeunits(channel),1),
            Ref(codeunits(""),1),
            1.0, 1000.0,
            units == "seconds" ? Val_Seconds : Val_Ticks,
            edge1 == "rising" ? Val_Rising : Val_Falling,
            edge2 == "rising" ? Val_Rising : Val_Falling,
            Ref(codeunits(""),1)) )
    t
end

"""
`generate_pulses(channel; low=2,high=2,delay=0) -> task`

create a NIDAQ counter output channel
"""
function generate_pulses(channel::String; low::T=2, high::T=2, delay::T=0) where T<:Number
    t = COTask()
    if T<:AbstractFloat
        ret = CreateCOPulseChanTime(t.th,
                Ref(codeunits(channel),1),
                Ref(codeunits(""),1),
                Val_Seconds,
                Val_Low,
                convert(Float64,delay),
                convert(Float64,low),
                convert(Float64,high))
    elseif T<:Integer
        ret = CreateCOPulseChanTicks(t.th,
                Ref(codeunits(channel),1),
                Ref(codeunits(""),1),
                Ref(codeunits(""),1),
                Val_Low,
                convert(Int32,delay),
                convert(Int32,low),
                convert(Int32,high))
    else
        error("low, high, and delay must either be \"FloatingPoint\" or \"Integer\"")
    end
    catch_error(ret)
    t
end

"""`read(task,channel; nsamples=-1) -> Vector[, Vector]`

receive data from the specified counter channel in the specified NIDAQ task
"""
function read(t::CITask, channel::String; num_samples::Integer = -1)
    if num_samples==-1;  num_samples=1024;  end

    #function read_counter_scalar(precision::DataType, cfunction::Function)
    #    data = precision[0]
    #    ret = cfunction(t, 1.0, pointer(data), pointer(C_NULL))
    #    ret>0 && @warn("NIDAQmx: $ret")
    #    ret<0 && error("NIDAQmx: $ret")
    #    data
    #end

    function read_counter_vector(precision::DataType, cfunction::Function)
        num_samples_read = Int32[0]
        data = Array(precision, num_samples)
        catch_error( cfunction(t.th,
            convert(Int32,num_samples),
            1.0,
            Ref(data,1),
            convert(UInt32,num_samples),
            Ref(num_samples_read,1),
            reinterpret(Ptr{Bool32},C_NULL)) )
        data = data[1:num_samples_read[1]]
        reshape(data, (num_samples, div(length(data),num_samples)))
    end

    function read_counter_2vectors(precision::DataType, cfunction::Function)
        num_samples_read = Int32[0]
        high = Array(precision, num_samples)
        low = Array(precision, num_samples)
        catch_error( cfunction(t.th,
            convert(Int32,num_samples),
            1.0,
            Val_GroupByChannel,
            Ref(high,1),
            Ref(low,1),
            convert(UInt32,num_samples),
            Ref(num_samples_read,1),
            reinterpret(Ptr{Bool32},C_NULL)) )
        high = high[1:num_samples_read[1]]
        low = low[1:num_samples_read[1]]
        reshape(high, (num_samples, div(length(high),num_samples))),
            reshape(low, (num_samples, div(length(low),num_samples)))
    end

    tmp = channel_type(t.th, channel)
    if tmp[2] == Val_CountEdges
        data = read_counter_vector(UInt32, ReadCounterU32)
    elseif tmp[2] == Val_PulseTime
        data = read_counter_2vectors(Float64, ReadCtrTime)
    elseif tmp[2] == Val_PulseTicks
        data = read_counter_2vectors(UInt32, ReadCtrTicks)
    elseif tmp[2] == Val_Position_AngEncoder
        data = read_counter_vector(UInt32, ReadCounterU32)
    elseif tmp[2] == Val_TwoEdgeSep  # might be broken
        val = Cint[0]
        catch_error( GetCITwoEdgeSepUnits(t.th,
                Ref(codeunits(channel),1),
                Ref(val,1)) )
        if val[1] == Val_Ticks
            data = read_counter_vector(UInt32, ReadCounterU32)
            #data = read_counter_scalar(UInt32, ReadCounterScalarU32)
        elseif val[1] == Val_Seconds
            data = read_counter_vector(Float64, ReadCounterF64)
            #data = read_counter_scalar(Float64, ReadCounterScalarF64)
        end
    end
    data
end


text:
using VarStructs: Props

    # Call Syntax
    person = @var Person(
                name::String = "Amin",
                number::Float64 = 20.0,
            )

    animal = @var Animal(
                name = "lion",
                number = 20,
            )

    # getproperty
    @test person.name == "Amin"
    @test person.number == 20.0

    # setproperty
    person.name = "Tim"
    person.number = 21.0
    @test person.name == "Tim"
    @test person.number == 21.0

    # fieldnames
    @test fieldnames(person) == [:name, :number]

    # add fields
    person.initial = "TH"
    @test fieldnames(person) == [:initial, :name, :number]

    # Dispatch
    function info(x::Person)
        println("Their home is city")
        return x.name
    end

    function info(x::Animal)
        println("Their home is jungle")
        return x.name
    end

    @test info(person) == "Tim"
    @test info(animal) == "lion"

    # getting an instance
    person = Person(name = "Amin", number = 20.0)

    # alternative way of getting instance
    person2 = Person(Dict(
        :name => Props("Amin"),
        :number => Props(20.0),
    ))

    @test info(person) == "Amin"
    @test info(animal) == "lion"

    # custom constructor
    function Person(name, number)
        return Person(Dict(
            :name => Props(name),
            :number => Props(number),
            :initial => Props(name[1]),
        ))
    end

    Person("Amin", 20.0)

    # Struct Syntax

    person2 = @var struct Person2
            name::String
            number::Float64 = 20.0
        end

    animal2 = @var struct Animal2
                name::String
                number::Int64
            end

    # getproperty
    @test isunset(person2.name)
    @test person2.number == 20.0

    # setproperty
    person2.name = "Tim"
    person2.number = 21.0
    @test person2.name == "Tim"
    @test person2.number == 21.0

    # fieldnames
    @test fieldnames(person2) == [:name, :number]

    # add fields
    person2.initial = "TH"
    @test fieldnames(person2) == [:initial, :name, :number]

    # Dispatch
    function info(x::Person2)
        println("Their home is city")
        return x.name
    end

    function info(x::Animal2)
        println("Their home is jungle")
        return x.name
    end

    # getting an instance
    person2 = Person2(name = "Amin", number = 20.0)
    animal2 = Animal2(name = "lion", number = 2)

    # alternative way of getting instance
    person3 = Person2(Dict(
        :name => Props("Amin"),
        :number => Props(20.0),
    ))

    @test info(person2) == "Amin"
    @test info(animal2) == "lion"

    mystruct = @var struct MyStruct
        x::Int64
    end

    @test isunset(mystruct.x)


text:
using QuantumStateDB
using Documenter

DocMeta.setdocmeta!(QuantumStateDB, :DocTestSetup, :(using QuantumStateDB); recursive=true)

makedocs(;
    modules=[QuantumStateDB],
    authors="JingYu Ning <foldfelis@gmail.com> and contributors",
    repo="https://github.com/foldfelis-QO/QuantumStateDB.jl/blob/{commit}{path}#{line}",
    sitename="QuantumStateDB.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://foldfelis-QO.github.io/QuantumStateDB.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/foldfelis-QO/QuantumStateDB.jl",
    devbranch="master",
)


text:
using JuLIPMaterials, ForwardDiff
using Test, JuLIP, JuLIP.Potentials, Einsum
using JuLIPMaterials.CLE, JuLIPMaterials.Testing
using JuLIPMaterials: Vec3, Mat3, ForceConstantMatrix1
using GaussQuadrature: legendre
using JuLIPMaterials.CLE: _C2, _dC2, _C2inv, _dC2inv, _dC2inv1, _ddC2inv,
                          _ddC2inv1, _C4, _dC4, _ddC4, _corrector_multiplier,
                          _d_corrector_multiplier, _dd_corrector_multiplier,
                          elastic_moduli

const CLE = JuLIPMaterials.CLE

using LinearAlgebra, Printf
##

# Set up a simple reference atomistic calculator
at = bulk(:Cu)
r0 = rnn(:Cu)
lj = lennardjones(r0=r0, rcut=[1.3*r0, 1.7*r0])
set_calculator!(at, lj)
variablecell!(at)
minimise!(at)
# Get force constants
fcm = ForceConstantMatrix1(lj, at, h = 1e-5)
# Get moduli
ℂ = elastic_moduli(at)

println("----------------------------------------------------------")
println(" Testing the 3D Green's Function Corrector Implementation")
println("----------------------------------------------------------")

@info("Test multiplier functions: dC2")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇C2 = ForwardDiff.derivative(t->_C2(x+t*a,ℂ),0.0)
    dC2 = _dC2(x,ℂ)
    dC2a = zeros(3,3)
    @einsum dC2a[i,j] = dC2[i,j,k]*a[k]
    global maxerr = max(maxerr, norm(dC2a-∇C2));
end
println("maxerr = $maxerr: ")
println(@test maxerr < 1e-9)

@info("Test multiplier functions: dC2inv")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇C2inv = ForwardDiff.derivative(t->_C2inv(x+t*a,ℂ),0.0)
    dC2inv = _dC2inv(x,ℂ)
    dC2inva = zeros(3,3)
    @einsum dC2inva[i,j] = dC2inv[i,j,k]*a[k]
    global maxerr = max(maxerr, norm(dC2inva-∇C2inv));
end
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-9)

@info("Test multiplier functions: dC2inv1")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇C2inv = ForwardDiff.derivative(t->_C2inv(x+t*a,ℂ),0.0)
    dC2inv1 = _dC2inv1(x,ℂ)
    dC2inv1a = zeros(3,3)
    @einsum dC2inv1a[i,j] = dC2inv1[i,j,k]*a[k]
    global maxerr = max(maxerr, norm(dC2inv1a-∇C2inv));
end
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-9)

@info("Test multiplier functions: ddC2inv")
@warn("THIS TEST IS FAILING; TODO -> FIX IT")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇²C2inv = ForwardDiff.derivative(s -> ForwardDiff.derivative(t->_C2inv(x+t*a,ℂ),s),0.0)
    ddC2inv = _ddC2inv(x,ℂ)
    ddC2invaa = zeros(3,3)
    @einsum ddC2invaa[i,j] = ddC2inv[i,j,k,l]*a[k]*a[l]
    global maxerr = max(maxerr, norm(ddC2invaa-∇²C2inv))
end
print("maxerr = $maxerr: ")
# println(@test maxerr < 1e-5)

@info("Test multiplier functions: ddC2inv1")
@warn("THIS TEST IS FAILING; TODO -> FIX IT")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇²C2inv = ForwardDiff.derivative(s -> ForwardDiff.derivative(t->_C2inv(x+t*a,ℂ),s),0.0)
    ddC2inv1 = _ddC2inv1(x,ℂ)
    ddC2inv1aa = zeros(3,3)
    @einsum ddC2inv1aa[i,j] = ddC2inv1[i,j,k,l]*a[k]*a[l]
    global maxerr = max(maxerr, norm(ddC2inv1aa-∇²C2inv));
end
print("maxerr = $maxerr: ")
# println(@test maxerr < 1e-5)

@info("Test multiplier functions: dC4")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇C4 = ForwardDiff.derivative(t->_C4(x+t*a,fcm),0.0)
    dC4 = _dC4(x,fcm)
    dC4a = zeros(3,3)
    @einsum dC4a[i,j] = dC4[i,j,k]*a[k]
    global maxerr = max(maxerr, norm(dC4a-∇C4));
end
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-9)

@info("Test multiplier functions: ddC4")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇²C4 = ForwardDiff.derivative(s -> ForwardDiff.derivative(t->_C4(x+t*a,fcm),s),0.0)
    ddC4 = _ddC4(x,fcm)
    ddC4aa = zeros(3,3)
    @einsum ddC4aa[i,j] = ddC4[i,j,k,l]*a[k]*a[l]
    global maxerr = max(maxerr, norm(ddC4aa-∇²C4));
end
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-5)

@info("Test multiplier functions: d_corrector_multiplier")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇H4 = ForwardDiff.derivative(t->_corrector_multiplier(x+t*a,ℂ,fcm),0.0)
    dH4 = _d_corrector_multiplier(x,ℂ,fcm)
    dH4a = zeros(3,3)
    @einsum dH4a[i,j] = dH4[i,j,k]*a[k]
    global maxerr = max(maxerr, norm(dH4a-∇H4));
end
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-9)

@info("Test multiplier functions: dd_corrector_multiplier")
@warn("TODO: THIS TEST IS FAILING")
maxerr = 0.0
for i=1:10
    a, x = Vec3(randvec3()), Vec3(randvec3())
    ∇²H4 = ForwardDiff.derivative(s->ForwardDiff.derivative(t->_corrector_multiplier(x+t*a,ℂ,fcm),s),0.0)
    ddH4 = _dd_corrector_multiplier(x,ℂ,fcm)
    ddH4aa = zeros(3,3)
    @einsum ddH4aa[i,j] = ddH4[i,j,k,l]*a[k]*a[l]
    global maxerr = max(maxerr, norm(ddH4aa-∇²H4));
end
println("maxerr = $maxerr")
# @test maxerr < 1e-5

# Set up Green's function and corrector
G0 = GreenFunction(ℂ, Nquad = 32)
Gcorr = GreenFunctionCorrector(ℂ, fcm, Nquad = 32)

println("Test that Gcorr satisfies the PDE div(C2 ∇Gcorr) = C4[∇]G0 ")
@info("TODO: THIS TEST IS FAILING!")
maxerr = 0.0
for i=1:5
    a, x = randvec3(), randvec3()
    # Compute LHS via AD
    v = x -> Gcorr(x)*a
    LHS = cleforce(Vec3(x),v,ℂ)
    # Compute RHS via AD
    RHS = zeros(3)
    for i=1:length(fcm.R)
        # Function to be differentiated
        u0 = t -> G0(x + t * fcm.R[i])*a
        # Compute 4th order derivative in ρ direction
        u1 = t -> ForwardDiff.derivative(u0,t)
        u2 = t -> ForwardDiff.derivative(u1,t)
        u3 = t -> ForwardDiff.derivative(u2,t)
        u4 = t -> ForwardDiff.derivative(u3,t)
        RHS += fcm.H[i]*u4(0.0)/24
    end
    # Compare
    global maxerr = max( norm(LHS-RHS, Inf), maxerr )
end
print("maxerr = $maxerr: ")
# println(@test maxerr < 1e-4)

println("Convergence of Gcorr (please test this visually!): ")
println(" nquad |    err    ")
println("-------|-----------")
xtest = [ randvec3() for n = 1:10 ]
Gref = GreenFunctionCorrector(ℂ, fcm, Nquad = 64)
for nquad in [2, 4, 8, 16, 32]
   global Gref, xtest, fcm, ℂ
   G = GreenFunctionCorrector(ℂ, fcm, Nquad = nquad)
   err = maximum( norm(Gref(x) - G(x), Inf)   for x in xtest )
   @printf("   %2d  | %.3e \n", nquad, err)
end


text:
module TriangleMeshes
using StaticArrays
using FastClosures
using DataStructures # heappush!,
using IGLWrap_jll
# libiglwrap="/home/jerome/src/iglwrap/local/libiglwrap.so"
# the bits types are hard-coded on the C side:
const Point=SVector{3,Cdouble}
const Face=NTuple{3,Cint}
@inline _face(a)=(Cint.(a[1:3])...,)

const Vec3d = MVector{3,Cdouble}
const Mat3d = MMatrix{3,3,Cdouble,9}

include("CornerTables.jl")
using .CornerTables

# Geometry««1
@inline det(a::SVector{2}, b::SVector{2}) = a[1]*b[2]-a[2]*b[1]
@inline det(a::SVector{2}, b::SVector{2}, c::SVector{2}) = det(a-c, b-c)
@inline distance2(a::SVector{3}, b::SVector{3}) = norm2(a-b)
@inline norm2(a::SVector{3}) = a[1]^2+a[2]^2+a[3]^2
# Data type ««1

struct TriangleMesh{T,A}
	vertices::Vector{SVector{3,T}}
	faces::Vector{Face}
	attributes::Vector{A}
	@inline TriangleMesh{T,A}(v, f, a) where{T,A} =
		new{T,A}(v, _face.(f), a)
	@inline TriangleMesh(v::AbstractVector{<:AbstractVector{T}}, f,
		a::AbstractVector{A}) where{T,A} = TriangleMesh{T,A}(v, f, a)
# 	@inline TriangleMesh{T}(v, f, a::AbstractVector{A}) where{T,A} =
# 		TriangleMesh{T,A}(v, f, a)
end

const CTriangleMesh = TriangleMesh{Cdouble}

@inline vertices(m::TriangleMesh) = m.vertices
@inline faces(m::TriangleMesh) = m.faces
@inline attributes(m::TriangleMesh) = m.attributes
@inline nvertices(m::TriangleMesh) = size(m.vertices, 1)
@inline nfaces(m::TriangleMesh) = size(m.faces, 1)
@inline shift(f::Face, k) = f .+ Face((k,k,k))

# IGL interface ««1

# a type to capture what is returned by IGL functions ««
@inline unsafe_array(T, p, n) =
	unsafe_wrap(Array, convert(Ptr{T}, p), n; own=true)
mutable struct Cmesh
	nv::Cint
	nf::Cint
	vertices::Ptr{Cdouble}
	faces::Ptr{Cint}
	@inline Cmesh() = new()
end
@inline fieldpointer(m::T, name) where{T} =
	pointer_from_objref(m) + fieldoffset(T, findfirst(==(name), fieldnames(T)))
@inline nvertices(m::Cmesh) = fieldpointer(m, :nv)
@inline nfaces(m::Cmesh) = fieldpointer(m, :nf)
@inline vpointer(m::Cmesh) = fieldpointer(m, :vertices)
@inline fpointer(m::Cmesh) = fieldpointer(m, :faces)
CTriangleMesh(m::Cmesh, a::AbstractVector{A}) where{A} = CTriangleMesh{A}(
	unsafe_array(Point, m.vertices, m.nv),
	unsafe_array(Face , m.faces   , m.nf),
	a)

@inline vpointer(m::TriangleMesh{Cdouble}) =
	convert(Ptr{Cdouble}, pointer(m.vertices))
@inline fpointer(m::TriangleMesh) = convert(Ptr{Cint}, pointer(m.faces))

#»»

function boolean(op, m1::CTriangleMesh{A}, m2::CTriangleMesh{A}) where{A}#««
	n = nfaces(m1)
	cm = Cmesh()
	j = Ref(Ptr{Cint}(0));
	r = @ccall libiglwrap.mesh_boolean(op::Cint,
		nvertices(m1)::Cint, nfaces(m1)::Cint,
		vpointer(m1)::Ref{Cdouble}, fpointer(m1)::Ref{Cint},
		nvertices(m2)::Cint, nfaces(m2)::Cint,
		vpointer(m2)::Ref{Cdouble}, fpointer(m2)::Ref{Cint},
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}},
		j::Ref{Ptr{Cint}})::Cint
	index = unsafe_array(Cint, j[], cm.nf)
	ao = [ i ≤ n ? m1.attributes[i] : m2.attributes[i-n] for i in index ]
	return CTriangleMesh(cm, ao)
end#»»
function minkowski_sum(m1::CTriangleMesh{A}, m2::CTriangleMesh{A}) where{A}#««
	n = nfaces(m1)
	cm = Cmesh()
	j = Ref(Ptr{Cint}(0));
	r = @ccall libiglwrap.minkowski_sum(
		nvertices(m1)::Cint, nfaces(m1)::Cint,
		vpointer(m1)::Ref{Cdouble}, fpointer(m1)::Ref{Cint},
		nvertices(m2)::Cint, nfaces(m2)::Cint,
		vpointer(m2)::Ref{Cdouble}, fpointer(m2)::Ref{Cint},
		3::Cint,
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}},
		j::Ref{Ptr{Cint}})::Cint
# 	index = unsafe_array(Cint, j[], cm.nf)
# 	ao = fill(first(attributes(m1)), size(index))
	return CTriangleMesh(cm, fill(first(m1.attributes), cm.nf))
# 	ao = [ i ≤ n ? m1.attributes[i] : m2.attributes[i-n] for i in index ]
# 	return (CTriangleMesh{A}(rvo, rfo, ao), index)
end#»»
function minkowski_sum(m1::CTriangleMesh{A}, v2::Vector{Point},
		e2::Vector{NTuple{2,Cint}}) where{A}#««
	n = nfaces(m1)
	cm = Cmesh()
	j = Ref(Ptr{Cint}(0));
	r = @ccall libiglwrap.minkowski_sum(
		nvertices(m1)::Cint, nfaces(m1)::Cint,
		vpointer(m1)::Ref{Cdouble}, fpointer(m1)::Ref{Cint},
		length(v2)::Cint, length(e2)::Cint,
		pointer(v2)::Ptr{Cdouble}, pointer(e2)::Ptr{NTuple{2,Cint}},
		2::Cint,
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}},
		j::Ref{Ptr{Cint}})::Cint
	# index = unsafe_array(Cint, j[], cm.nf)
# 	ao = fill(first(attributes(m1)), size(index))
# 	return (rvo, rfo, index)
	return CTriangleMesh(cm, fill(first(attributes(m1)), cm.nf))
# 	ao = [ i ≤ n ? m1.attributes[i] : m2.attributes[i-n] for i in index ]
# 	return (CTriangleMesh{A}(rvo, rfo, ao), index)
end#»»
function ispwn(m::TriangleMesh)#««
	r = ccall((:mesh_is_pwn, libiglwrap), Cint,
		(Cint, Cint, Ref{Cdouble}, Ref{Cint},),
		nvertices(m), nfaces(m), vpointer(m), fpointer(m),
		)
	return (r ≠ 0)
end#»»
function offset(m::CTriangleMesh{A}, level::Real, grid::Integer) where{A}#««
	cm = Cmesh()
	r = @ccall libiglwrap.offset_surface(
		nvertices(m)::Cint, nfaces(m)::Cint,
		vpointer(m)::Ref{Cdouble}, fpointer(m)::Ref{Cint},
		level::Cdouble, grid::Cint,
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}})::Cint
	@assert r == 0

	return CTriangleMesh(cm, fill(first(m.attributes), cm.nf))
end#»»
function decimate(m::CTriangleMesh{A}, max_faces::Integer) where{A}#««
	cm = Cmesh()
	j = Ref(Ptr{Cint}(0));
	r = @ccall libiglwrap.decimate(
		nvertices(m)::Cint, nfaces(m)::Cint,
		vpointer(m)::Ref{Cdouble}, fpointer(m)::Ref{Cint},
		max_faces::Cint,
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}},
		j::Ref{Ptr{Cint}})::Cint
	@assert r == 0

	index = unsafe_array(Cint, j[], cm.nf)
	return CTriangleMesh(cm, m.attributes[index])
end#»»
function loop(m::CTriangleMesh{A}, count::Integer) where{A}#««
	cm = Cmesh()
	r = @ccall libiglwrap.loop(
		nvertices(m)::Cint, nfaces(m)::Cint,
		vpointer(m)::Ref{Cdouble}, fpointer(m)::Ref{Cint},
		count::Cint,
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}})::Cint
	@assert r == 0

	return CTriangleMesh(cm, attributes(m)[fld1.(1:cm.nf, 4^count)])
end#»»
function centroid_and_volume(m::CTriangleMesh)#««
	c = Vec3d(undef)
	v = @ccall libiglwrap.centroid_and_volume(
		nvertices(m)::Cint, nfaces(m)::Cint,
		vpointer(m)::Ref{Cdouble}, fpointer(m)::Ref{Cint},
		c::Ref{Vec3d})::Cdouble
	return (c, v)
end#»»
function halfspace(direction, origin, m::CTriangleMesh{A}, color::A) where{A}#««
	cm = Cmesh()
	j = Ref(Ptr{Cint}(0));
	r = @ccall libiglwrap.intersect_with_half_space(
		nvertices(m)::Cint, nfaces(m)::Cint,
		vpointer(m)::Ref{Cdouble}, fpointer(m)::Ref{Cint},
		Vec3d(origin...)::Ref{Vec3d},
		Vec3d(direction...)::Ref{Vec3d},
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}},
		j::Ref{Ptr{Cint}})::Cint
	@assert r == 0
	index = unsafe_array(Cint, j[], cm.nf)

	return CTriangleMesh(cm,
		[(i <= nfaces(m) ? m.attributes[i] : color) for i in index])
end#»»
function swept_volume(m::CTriangleMesh, transform, #««
	nsteps, gridres, isolevel = 0)
	cm = Cmesh()
	f = (t,a,b) -> begin
		u, v = transform(t)
		unsafe_store!(a, Mat3d(u)) # column-major, same as Eigen
		unsafe_store!(b, Vec3d(v))
		return
	end
	ctransform = @cfunction($f, Cvoid, (Cdouble, Ptr{Mat3d}, Ptr{Vec3d}))
	r = @ccall libiglwrap.swept_volume(
		nvertices(m)::Cint, nfaces(m)::Cint,
		vpointer(m)::Ref{Cdouble}, fpointer(m)::Ref{Cint},
		ctransform::Ptr{Cvoid}, nsteps::Csize_t, gridres::Csize_t,
		isolevel::Cdouble,
		nvertices(cm)::Ptr{Cint}, nfaces(cm)::Ptr{Cint},
		vpointer(cm)::Ptr{Ptr{Cdouble}}, fpointer(cm)::Ptr{Ptr{Cint}})::Cint
	@assert r == 0

	return CTriangleMesh(cm, fill(first(attributes(m)), cm.nf))
end#»»

@inline Base.union(m1::CTriangleMesh, m2::CTriangleMesh) = boolean(0, m1, m2)
@inline Base.intersect(m1::CTriangleMesh, m2::CTriangleMesh)= boolean(1, m1, m2)
@inline Base.setdiff(m1::CTriangleMesh, m2::CTriangleMesh) = boolean(2, m1, m2)
@inline Base.symdiff(m1::CTriangleMesh, m2::CTriangleMesh) = boolean(3, m1, m2)

# Own functions ««1
# plane_splice ««2
"""
    plane_slice(m::TriangleMesh)

Returns the set of all edges formed by this mesh ∩ the horizontal plane,
as `(vertices, edges)`, where `vertices` are 2d points,
and `edges` are indices into `vertices`.
"""
function plane_slice(z::Real, m::TriangleMesh)
	# build a list of intersection points + connectivity
	# each intersection point is either:
	#  - a vertex v, represented as (v, 0)
	#  - in the edge (v1v2), represented as (v1,v2)
	points = Dict{NTuple{2,Int},Int}()
	elist = NTuple{2,Int}[]
	pindex = @closure v->get!(points, extrema(v), length(points)+1)
	edge! = @closure (v,w)-> push!(elist, minmax(pindex(v),pindex(w)))
	# build list of all edges:
	for (i1, i2, i3) in faces(m)
		(v1, v2, v3) = vertices(m)[[i1,i2,i3]]
		(z1, z2, z3) = (v1[3], v2[3], v3[3]) .- z
		if z1 == 0#««
			if z2 == 0
				z3 == 0 && continue # 000: ignore horizontal triangle
				edge!(i1=>0,i2=>0) # 00+, 00-
			elseif z2 > 0
				z3 == 0 && edge!(i1=>0,i3=>0)
				z3 < 0 && edge!(i1=>0,i2=>i3)
			else # z2 < 0
				z3 == 0 && edge!(i1=>0,i3=>0)
				z3 > 0 && edge!(i1=>0,i2=>i3)
			end
		elseif z1 > 0
			if z2 == 0
				z3 == 0 && edge!(i2=>0,i3=>0) # +00
				z3 < 0 && edge!(i2=>0,i1=>i3) # +0-
			elseif z2 > 0
				z3 < 0 && edge!(i1=>i3,i2=>i3) # ++-
			else # z2 < 0
				z3 == 0 && edge!(i3=>0,i1=>i2)
				z3 > 0 && edge!(i1=>i2,i2=>i3)#+-+
				z3 < 0 && edge!(i1=>i2,i1=>i3)#+--
			end
		else # z1 < 0
			if z2 == 0
				z3 == 0 && edge!(i2=>0,i3=>0) # -00
				z3 > 0 && edge!(i2=>0,i1=>i3) # -0+
			elseif z2 > 0
				z3 == 0 && edge!(i3=>0,i1=>i2)
				z3 > 0 && edge!(i1=>i2,i1=>i3)#-++
				z3 < 0 && edge!(i1=>i2,i2=>i3)#-+-
			else # z2 < 0
				z3 > 0 && edge!(i1=>i3,i2=>i3) # --+
			end
		end#»»
	end
	vlist = Vector{SVector{2,Float64}}(undef, length(points))
	for ((i1,i2),j) in pairs(points)
		if iszero(i1)
			v = vertices(m)[i2]
			vlist[j] = SA[v[1],v[2]]
		else
			(v1,v2) = vertices(m)[[i1,i2]]
			(z1,z2) = (v1[3]-z,v2[3]-z); f = 1/(z2-z1)
			# (z2 v1 - z1 v2)/(z2-z1)
			vlist[j] = SA[f*(z2*v1[1]-z1*v2[1]), f*(z2*v1[2]-z1*v2[2])]
		end
	end
	return (vlist, unique!(sort!(elist)))
end

# project««2
"""
    project(mesh)

Returns the list of triangles produced as the 2d projection of `mesh`.
"""
function project(m::TriangleMesh{T}) where{T}
	triangles = SVector{3,SVector{2,T}}[]
	for (i,j,k) in faces(m)
		a,b,c = vertices(m)[SA[i,j,k]]
		a1,b1,c1 = (SA[a[1],a[2]], SA[b[1],b[2]], SA[c[1],c[2]])
		d = det(a1,b1,c1)
		iszero(d) && continue
		push!(triangles, d > 0 ? SA[a1,b1,c1] : SA[a1,c1,b1])
	end
	return triangles
end

# Corner table structure ««2
struct CTMesh{J,T,A} <: AbstractTriangulation{J}
	triangulation::CornerTable{J}
	points::Vector{SVector{3,T}}
	attributes::Vector{A}
end

Base.show(io::IO, t::CTMesh) = CornerTables.showall(io, t)

@inline CornerTables.triangulation(t::CTMesh) = t.triangulation
@inline point(t::CTMesh, c::Cell) = t.points[int(c)]
@inline edgelength(t::CTMesh, a::Arrow; distance2 = distance2) =
	distance2(point(t, CornerTables.head(t,a)), point(t, CornerTables.tail(t,a)))

# conversion to and from `TriangleMesh`

@inline CTMesh{J}(m::TriangleMesh{T,A}) where{J,T,A} =
	CTMesh{J,T,A}(CornerTable{J}(faces(m)), copy(vertices(m)), copy(attributes(m)))
@inline CTMesh(m::TriangleMesh) = CTMesh{Int32}(m)

@inline TriangleMesh(t::CTMesh{J,T,A}) where{J,T,A} =
	TriangleMesh(t.points, [ Int.(f) for f in alltriangles(t) ], t.attributes)

# splitedges ««2
# The following structure is used as a heap for the edges of the mesh,
# sorted by size; we will always remove (and split) the longest edge:
struct VecSortedSet{J,V}
	set::SortedSet{Tuple{V,J}}
	size::Vector{V}
end
@inline VecSortedSet{J}(s::AbstractVector{V}) where{J,V} =
	VecSortedSet{J,V}(SortedSet{Tuple{V,J}}((x,i) for (i,x) in pairs(s)), s)
@inline VecSortedSet(s::AbstractVector) = VecSortedSet{Int}(s)

@inline Base.last(s::VecSortedSet) = last(s.set)[2]
@inline function Base.setindex!(s::VecSortedSet, v, k)
	delete!(s.set, (s.size[k], k))
	s.size[k] = v
	push!(s.set, (s.size[k], k))
	return s
end
@inline function Base.push!(s::VecSortedSet, vlist...)
	for v in vlist
		push!(s.size, v)
		push!(s.set, (last(s.size), length(s.size)))
	end
	return s
end
"""
    splitedges!(m::TriangleMesh, maxlen)

Repeatedly split all edges of `m`, starting by the longest ones,
until no edge has length² > `maxlen`.
"""
@inline splitedges(m::TriangleMesh, maxlen; distance2=distance2) =
	TriangleMesh(splitedges!(CTMesh(m), maxlen; distance2))

function splitedges!(t::CTMesh{J,T,A}, maxlen; distance2=distance2) where{J,T,A}
	elist = VecSortedSet{J}(T[ a > opposite(t, a) ? zero(T) :
		edgelength(t, a; distance2) for a in eacharrow(t) ])
	while true
	e = Arrow(last(elist)); o = opposite(t, e)
	elist.size[e] ≤ maxlen && break
	ne, po = next(e), prev(o)
	o_ne, o_po = opposite(t, ne), opposite(t, po)

	ct, ch = CornerTables.tail(t, e), CornerTables.head(t, e)
	pm = (point(t,ct) + point(t,ch))/2
	push!(t.points, pm)
	push!(t.attributes, t.attributes[int(node(e))], t.attributes[int(node(o))])
	(q1, q2, cm) = CornerTables.split!(t, e)
	# update elist accordingly:
	# edge (ct, cm) is now (e, opp(e))  (with e < opp(e))
	# edge (cm, ch) is (side(q1,1) < side(q2,1))
	# edge (cm, cl) is (next(e) < side(q1,3))
	# edge (cm, cr) is (prev(o) < side(q2,2))
	# edge (ch, cl) is (opp(next(e)) < side(q1,2))
	# edge (ch, cr) is (opp(prev(o)) < side(q2,3))
	pt, ph, pl, pr = point(t, ct), point(t, ch),
		point(t, left(t,e)), point(t, right(t,e))
	elist[e] = distance2(ph, pm) # elist[o] is already 0
	push!(elist, distance2(pm, ph), 0, 0, 0, 0, 0)
	elist[ne] = distance2(pm, pl)
	elist[po] = distance2(pm, pr)
	elist[o_ne] = distance2(ph, pl)
	elist[o_po] = distance2(ph, pr)
	end
	return t
end

#  »»1
export TriangleMesh
end


text:
## Some defaults if not already defined

DRT = isdefined(Main, :DRT) ? DRT : true
PPE = isdefined(Main, :PPE) ? PPE : true
REF = isdefined(Main, :REF) ? REF : true

##

dcx = 1200 # 1920
dcy = 771  # 1080

global scene
scene, layout = layoutscene(resolution = (dcx, dcy)) #10,

ax2 = LAxis(scene)
ax2.xlabel = "East [meters]"
ax2.ylabel = "North [meters]"

figaxes = [ax2]
tightlimits!.(figaxes)
layout[1:1, 1:1] = figaxes

##

scl =  1.0 #5.0
ofsx = 140.0    # 90*scl
ofsy = 170.0    # 145*scl

addLinesBelief!(fg, figaxes[1,1], TTm, scale=scl, origin=(ofsx,ofsy), drt=DRT, ppe=PPE, ref=REF, maskRef=(Second(420),Second(450)) )

scene
##

# resolution=(dcx,dcy)
# pl = plotVariableBeliefs(fg, r"x\d", sortVars=true, fade=2, fadeFloor=0.2, scene=figaxes[1,1], scale=scl, origin=(ofsx,750))
pl, Z = plotVariableBeliefs(fg, r"x\d", colormap=:blues, minColorBase=0.0,
                            sortVars=true, fade=1, extend=0.0,fadeFloor=0.2, scale=scl, origin=(ofsx,ofsy), N=100, scene=figaxes[1,1],
                            xmin=140-ofsx, xmax=240-ofsx,ymin=90-ofsy,ymax=130-ofsy, force = true)
                            # xmin=(0-ofsx)/scl, xmax=(dcx-ofsx)/scl,ymin=(0-ofsy)/scl,ymax=(dcy-ofsy)/scl)
#

scene

##


text:
# run this to make a precompiled system image with the packages used
# in Econometrics pre-compiled, and with some commonly used functions
# also pre-compiled.
#
# once it's done, start Julia from the Econometrics directory with
# julia --proj -JEconometricsSysimage.so"
#
# just start without the switch to return to normal. You will
# need to re-do this if you update your packages.

using Pkg, PackageCompiler
Pkg.activate(".")
PackageCompiler.create_sysimage(:Econometrics; sysimage_path="JuliaSysimage.so",
                                                precompile_execution_file="warmup.jl")

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
using Flight
using OrdinaryDiffEq
using SciMLBase
using LinearAlgebra
using BenchmarkTools

plot_settings = (linewidth=2, margin = 10mm, guidefontsize = 12)

trn = DummyTerrain()
atm_sys = System(AtmosphereDescriptor())

ac = AircraftBase(
    mass = TunableMass(),
    ctrl = NoControlMapping(),
    aero = NullAirframeComponent(),
    pwp = NullAirframeComponent(),
    ldg = NullAirframeComponent(),
    srf = NullAirframeComponent()
)

ac_sys = System(ac); #should remake the system, because it sets the Model's initial condition upon creation
init!(ac_sys.x.kin, KinInit(v_eOb_b = [10, 0, 0], Ob = Geographic(LatLon(ϕ = 0), AltOrth(1000))))
ac_mdl = Model(ac_sys, (trn, atm_sys); dt = 0.01, adaptive = false, solver = Heun(), y_saveat = 0.1);
solve!(ac_mdl)
plots(ac_mdl; save_path = joinpath("tmp", "plots_ac_01"), plot_settings...)

x0 = copy(ac_mdl.x)
init!(x0.kin, KinInit(v_eOb_b = [10, 0, 0], ω_lb_b = [0.5, 0.5, 0], Ob = Geographic(LatLon(ϕ = 0), AltOrth(1000))))
reinit!(ac_mdl, x0)
solve!(ac_mdl)
plots(ac_mdl; save_path = joinpath("tmp", "plots_ac_02"), plot_settings...)

################################## AC2 ######################################

ac = AircraftBase(
    kin = KinLTF(),
    mass = TunableMass(),
    aero = SimpleDrag(),
    pwp = AirframeGroup((
        left = EThruster(motor = ElectricMotor(α = CW)),
        right = EThruster(motor = ElectricMotor(α = CCW)))),
)

atm_sys.u.wind.v_ew_n[1] = -5
ac_sys = System(ac);
init!(ac_sys.x.kin, KinInit(v_eOb_b = [10, 0, 0], Ob = Geographic(LatLon(ϕ = 0), AltOrth(1000))))
ac_mdl = Model(ac_sys, (trn, atm_sys); dt = 0.01, adaptive = false, solver = Heun(), y_saveat = 0.1);
solve!(ac_mdl)
plots(ac_mdl; save_path = joinpath("tmp", "plots_ac_03"), plot_settings...)

ac_sys.subsystems.pwp.u.left.throttle = .1
ac_sys.subsystems.pwp.u.right.throttle = .1
reinit!(ac_mdl)
solve!(ac_mdl)
plots(ac_mdl; save_path = joinpath("tmp", "plots_ac_04"), plot_settings...)

ac_sys.subsystems.pwp.u.left.throttle = 0
ac_sys.subsystems.pwp.u.right.throttle = .2
reinit!(ac_mdl)
solve!(ac_mdl)
plots(ac_mdl; save_path = joinpath("tmp", "plots_ac_05"), plot_settings...)

b = @benchmarkable step!($ac_mdl, 1, true) setup=(reinit!($ac_mdl)); run(b)

# ac_sys = System(ac); #should remake the system, because it sets the Model's initial condition upon creation
# ac_mdl = Model(ac_sys, (trn, atm_sys));
# ac_sys.subsystems.pwp.u.left.throttle = 1 #the same
# b = @benchmarkable step!($ac_mdl, 1, true) setup=(reinit!($ac_mdl)); run(b)

# #kinematics torture test: do not try this without an adaptive integrator! the
# #"thing" shoots southward from the North pole. if you're looking exactly south,
# #your heading is ill-defined, so it chatters between π and -π at different
# #integration step. but the quaternion representation doesn't care. bank angle
# #wraps around -π twice longitude is not defined at the North pole, so it jumps
# #instantaneously from its random initialization value (0) to around π. it
# #doesn't chatter because due to Coriolis acceleration, it does not remain at π.

# #reinitialize at the North Pole, set end time and run to completion
# ac_sys.subsystems.pwp.u.left.throttle = 0. #the same
# ac_sys.subsystems.pwp.u.right.throttle = 0. #the same
# atm_sys.u.wind.v_ew_n[1] = -5
# reinit!(ac_mdl)
# x0 = copy(ac_mdl.x)
# init!(x0.kin, KinInit(v_eOb_b = [10, 0, 0], Ob = Geographic(LatLon(ϕ = 0), AltOrth(1000))))
# reinit!(ac_mdl, x0, tf = 10)
# solve!(ac_mdl)

# #this should be set at startup.jl
# plot_settings = (linewidth=2, margin = 10mm, guidefontsize = 12)

# plots(ac_mdl; save_path = joinpath("tmp", "plots_drag"), plot_settings...)


text:
module repl

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Distributions
using Soss
using MeasureTheory
using StatsPlots
using Parquet
using DataFrames
using Chain
using UUIDs
using DataFrameMacros
using GLM

mode = "test"
# mode = "prod"

if mode == "test"
    n_simulations = 75
    n_applications_per_period = 50
elseif mode == "prod"
    n_simulations = 30
    n_applications_per_period = 1000
end

function generate_synthetic_data(n_applications_per_period)
    simulation_id = string(UUIDs.uuid4())

    loan_data_generator = @model age_var begin

        # Age distribution (mean, var) is an increasing function of time, min age of population is set to age_var / 4 for a given time period
        age ~ Distributions.TruncatedNormal(age_var, age_var, age_var / 4, 100)

        # Idiosyncratic individual risk is always zero mean, but the variance is an inverse function of age
        # Older individuals are as such modeled as riskier, but more consistent; younger are safer, but less predictable
        idiosyncratic_individual_risk ~ Distributions.TruncatedNormal(0, 1 / (age * 5), -10, 10)

        # Age is assumed to be associated with increased risk; average risk is located at logistic(-3), approx 4.7% default risk, plus age effects
        total_default_risk_log_odds = idiosyncratic_individual_risk + age - 3

        # Log odds scale is mapped onto zero to one score
        total_default_risk = logistic(total_default_risk_log_odds)

        # Default is a weighted coin flip with weight of 'total_default_risk'
        default ~ MeasureTheory.Bernoulli(total_default_risk)
    end

    # Age varible is an increasing function of time, as specified here for each time step
    age_var = [0.25, 0.25, 0.3, 0.4, 0.5, 1, 1.2, 1.5, 2.0, 2.0]

    business_cycle_df = DataFrame(
        "age_var" => age_var,
    )

    n_periods = nrow(business_cycle_df)
    business_cycle_df[!, "application_date"] = 2020:(2020 + (n_periods - 1))

    # A separate sample is drawn for each moment in time, based on the age_var for that cross-section
    portfolio_df = DataFrame()
    for x in eachrow(business_cycle_df)
        one_cycle_df = DataFrame(rand(loan_data_generator(x.age_var), n_applications_per_period))
        one_cycle_df = @chain one_cycle_df begin
            @transform(:application_date = x.application_date,
                       :default = :default * 1, # convert bool to int
                       :age_var = x.age_var,
                       )
        end
        append!(portfolio_df, one_cycle_df)
    end

    # UUIDs are generated for easier tracking
    portfolio_df = @chain portfolio_df begin
        @transform(
            :application_id = string(UUIDs.uuid4()),
            :simulation_id = simulation_id,
        )
    end

    file_path = "/app/$(simulation_id).parquet"

    # Data is saved in parquet format
    @chain portfolio_df begin
        write_parquet(file_path, _)
    end

    # Data is sent to cloud
    bucket_name = ENV["S3_BUCKET_NAME"]
    run(`aws s3api put-object --bucket $(bucket_name) --key synthetic_data/$(simulation_id).parquet --body $(simulation_id).parquet`)
    rm(file_path)

    return portfolio_df
end

function generate_synthetic_data(n_applications_per_period, n_simulations)
    bucket_name = ENV["S3_BUCKET_NAME"]

    # Clear bucket of synthetic data...
    run(`aws s3 rm s3://$(bucket_name)/synthetic_data --recursive`)

    for i = 1:n_simulations
        generate_synthetic_data(n_applications_per_period)
    end
end

generate_synthetic_data(n_applications_per_period, n_simulations)

max_stars_repo_name: jeremiahpslewis/alec

max_issues_repo_name: jeremiahpslewis/alec

max_forks_repo_name: jeremiahpslewis/alec


text:
@testgroup "options" begin

@testcase "empty options" begin
    opts = Jute.build_run_options()

    # test some default option
    @test opts[:verbosity] == 1
end

@testcase "commandline options" begin
    opts = Jute.build_run_options(; args=["--verbosity=2"])
    @test opts[:verbosity] == 2
end

@testcase "userdict options" begin
    opts = Jute.build_run_options(; options=Dict(:verbosity => 2))
    @test opts[:verbosity] == 2
end

@testcase "commandline overrides userdict options" begin
    opts = Jute.build_run_options(; options=Dict(:verbosity => 2), args=["--verbosity=0"])
    @test opts[:verbosity] == 0
end

@testcase "commandline option name change" begin

    # Confirm the default option
    opts = Jute.build_run_options()
    @test !opts[:dont_add_runtests_path]

    # When the default value of a flag is set to true by userdict,
    # the flag is replaced by its inverse in the possible commandline options,
    # so that it could be disabled via commandline
    opts = Jute.build_run_options(;
        options=Dict(:dont_add_runtests_path => true),
        args=["--add-runtests-path"])
    @test !opts[:dont_add_runtests_path]

    opts = Jute.build_run_options(;
        options=Dict(:dont_add_runtests_path => false),
        args=["--dont-add-runtests-path"])
    @test opts[:dont_add_runtests_path]

end

@testcase "unknown option in userdict" begin
    @test_throws ErrorException Jute.build_run_options(; options=Dict(:foobar => 0))
end

@testcase "tag options processing" begin
    opts = Jute.build_run_options(; options=Dict(:include_only_tags => [:tag1, :tag2]))
    @test opts[:include_only_tags] == [:tag1, :tag2]

    opts = Jute.build_run_options(; args=["-t", "tag1", "tag2"])
    @test opts[:include_only_tags] == [:tag1, :tag2]

    opts = Jute.build_run_options(; options=Dict(:exclude_tags => [:tag1, :tag2]))
    @test opts[:exclude_tags] == [:tag1, :tag2]

    opts = Jute.build_run_options(; args=["-n", "tag1", "tag2"])
    @test opts[:exclude_tags] == [:tag1, :tag2]
end

@testcase "color options processing" begin
    opts = Jute.build_run_options(; options=Dict(:color_pass => :blue))
    @test opts[:report_color_scheme].color_pass == :blue

    opts = Jute.build_run_options(; args=["--color-pass=blue"])
    @test opts[:report_color_scheme].color_pass == :blue
end

end


text:
"""
    order_tw_by_dissections_threads(graph, duration,nparallel=10, max_imbalances=[1.0, 0.8, 0.6], seed=4242)
Computes an approximation of the upper bound of the treewidth of the graph g and
an order for elimination using the iterative dissection (using flow cutter algorithm)
several times in parallel and select the best_order and treewidth.

# Arguments
- `graph::SimpleGraph{Int64}`: the graph to analyse.
- `duration::Int64` : the maximum time to wait (in seconds) for running iterative_dissection multiple times.
- `nparallel::Int64 = 10` the number of calls in parallel to iterative_dissection.
- `max_imbalance::Float64 = 0.6` criteria for selecting cuts that will build the separator.
- `max_nsample::Int64 = 20` the number of calls to flowcutter with random inputs.
- `seed::Int64 = 4242` the base seed for the RNG sampling the inputs of flowcutter and the choosen one cut for separator.

# Return
- `order::Vector{Int64}` an array of vertices index.
- `treewidth::Int64` the approximation of treewidth.

"""
function order_tw_by_dissections_threads(
    graph::SimpleGraph{Int64},
    duration::Int64,
    nparallel::Int64 = 2 * Threads.nthreads(),
    max_imbalances::Vector{Float64} = [1.0, 0.8, 0.6],
    seed::Int64 = 4242
    )::Pair{Vector{Int64},Int64}

    best_tw = typemax(Int64)
    best_order = Vector{Int64}()
    start = time()
    k = length(max_imbalances)

    while time() - start < duration
        @debug "start $nparallel run current tw = $best_tw , time = $(time() - start)"
        seed += Threads.nthreads()
        Threads.@threads for i = 1:nparallel
            res = iterative_dissection(graph, best_tw, max_imbalances[i % k + 1] , 20, seed + i)
            if res[2] < best_tw
                best_order, best_tw = res
            end
        end
    end
    @debug "finally best tw = $best_tw"
    return best_order => best_tw
end

"""
    order_tw_by_dissections_simple(graph, duration,)
Computes an approximation of the upper bound of the treewidth of the graph g and
an order for elimination using the iterative dissection (using flow cutter algorithm)
several times and select the best_order and treewidth.

# Arguments
- `graph::SimpleGraph{Int64}`: the graph to analyse.
- `duration::Int64` : the maximum time to wait (in seconds) for running iterative_dissection multiple times.
- `max_imbalance::Float64 = 0.6` criteria for selecting cuts that will build the separator.
- `seed::Int64 = 4242` the base seed for the RNG sampling the inputs of flowcutter and the choosen one cut for separator.
- `nsample::Int64 = 20` the number of calls to flowcutter with random inputs.
- `sample_augment::Int64 = 8` the number of iteration required before augmenting nsample.

# Return
- `order::Vector{Int64}` an array of vertices index.
- `treewidth::Int64` the approximation of treewidth.

"""
function order_tw_by_dissections_simple(
    graph::SimpleGraph{Int64},
    duration::Int64;
    max_imbalances::Vector{Float64} = [1.0, 0.8, 0.6],
    seed::Int64 = 4242,
    nsample::Int64 = 20,
    sample_augment::Int64 = 8
    )::Pair{Vector{Int64},Int64}

    best_tw = typemax(Int64)
    best_order = Vector{Int64}()
    start = time()
    k = length(max_imbalances)
    i = 1

    while time() - start < duration
        @debug "time = $(round(time() - start; digits=2))"
        @debug "start new run with seed = $seed , current tw = $best_tw , max_imbalance = $(max_imbalances[i % k + 1]), nsample = $(20 + 20*(i÷8))"
        res = iterative_dissection(graph, best_tw, max_imbalances[i % k + 1] , nsample + nsample*(i ÷ sample_augment), seed)
        if res[2] < best_tw
            best_order, best_tw = res
        end
        seed += 1
        i += 1
    end
    @debug "finally best tw = $best_tw"
    return best_order => best_tw
end


text:
using DataStructures

"""
Hybrid Actor-critic
-------------------

Hybrid model : Actor-critic & Forward learner

See learner/forwardlearner.jl, learner/actorcriticlearner.jl
and learner/hybridlearner.jl for the core implementation.

Note 1
------
get_all_puzzles(), get_all_episodes_in_puzzle(), get_data_by_puzzle_episode()
are in rl_data_toolsforalgos.jl

Note 2
------
Temperature of softmax is removed, because this learner in NOT convex and weights
can be in (0, +Inf). The temperature is absorbed in the weights. Otherwise
it is overparametrized.
"""
function getLL_Hybrid_AC_FWD(data::Array{Int,2},
                            alpha::Float64, eta_critic::Float64,
                            gamma::Float64, lambda_actor::Float64,
                            lambda_critic::Float64,
                            eta_fwd::Float64,
                            w_MB::Float64,
                            w_MF::Float64;
                            init_R = 0.,
                            init_policy=nothing,
                            init_Vs = nothing,
                            use_replacing_trace=true,
                            do_debug_print = false)

    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_Hybrid_AC_FWD supports only one subject")
    end
    all_puzzles = get_all_puzzles(data)

    LL = 0.
    action_selection_probs = zeros(Float64, nr_samples)
    # signals:
    TD_record = zeros(Float64, nr_samples)
    SPE_record = zeros(Float64, nr_samples)
    V_state = zeros(Float64, nr_samples + data_stats.nr_rewards)

    PolicyParametersDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)
    PolicyDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)

    learnerMB = []
    learnerMF = []
    learnerHybrid = []
    modulatingfactor = 0.
    action_counter = 0
    for i_puzzle in all_puzzles
        # Puzzle = independent block of the experiment.
        # The code was made for general purposes, in this experiment we
        # have only ONE puzzle.
        learnerMB = ForwardLearner(nr_states = nr_states,
                            nr_actions = nr_actions,
                            η = eta_fwd, γ = gamma,
                            init_R = 0.
                            )
        if init_policy == nothing # providing init_policy w/o init_Vs is considered a bug
            learnerMF = ActorCritic(nr_states = nr_states,
                                nr_actions = nr_actions,
                                α_critic = eta_critic,
                                α_actor = alpha,
                                γ = gamma,
                                λ_critic = lambda_critic,
                                λ_actor = lambda_actor,
                                use_replacing_trace = use_replacing_trace,
                                do_debug_print = do_debug_print)
        else
            learnerMF = ActorCritic(nr_states = nr_states,
                                nr_actions = nr_actions,
                                α_critic = eta_critic,
                                α_actor = alpha,
                                γ = gamma,
                                λ_critic = lambda_critic,
                                λ_actor = lambda_actor,
                                policyparameters = copy(init_policy[current_puzzle_id]),
                                V = copy(init_Vs[current_puzzle_id]),
                                use_replacing_trace = use_replacing_trace,
                                do_debug_print = do_debug_print)
        end
        learnerHybrid = Hybrid(nr_states = nr_states,
                            nr_actions = nr_actions,
                            w_MB = w_MB,
                            w_MF = w_MF
                            )
        all_epi_in_puzzle = get_all_episodes_in_puzzle(data, i_puzzle)
        j = 1
        for i_episode in all_epi_in_puzzle
            episode_data = get_data_by_puzzle_episode(data, i_puzzle, i_episode)
            path_length = size(episode_data)[1]  # = T
            if do_debug_print
                println("i_puzzle:$i_puzzle | i_episode:$i_episode | path_length:$path_length")
            end
            isendofexperiment = false

            for t in 1:path_length
                action_counter += 1
                state = episode_data[t, Int(c_state)]
                action = episode_data[t, Int(c_action)]
                reward = episode_data[t, Int(c_reward)]
                next_state = episode_data[t, Int(c_nextState)]
                done = reward>0 ? true : false
                # @show state, action, reward, next_state
                if do_debug_print
                    # println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | p(s,a):$(round(p_action,2)) | p(s,1;+-bias): $(round(p1_pos,2))/$(round(p1_neg,2))")
                    println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | done: $done")
                end

                # Hybrid learner
                update!(learnerHybrid, state, action, reward, learnerMB.QSA, learnerMF.policyparameters) #Q_hybrid =  w_MB*learnerMB.QSA + w_MF*learnerMF.QSA

                p_actions = getSoftmaxProbs(learnerHybrid.QSA[state, :], 1., [0., 0.])
                action_selection_probs[action_counter] = copy(p_actions[action])
                LL += log(p_actions[action])

                PolicyParametersDiffChosenVsUnchosen_record[action_counter] = learnerHybrid.QSA[state, action] - learnerHybrid.QSA[state, mod(action,2)+1]
                PolicyDiffChosenVsUnchosen_record[action_counter] = p_actions[action] - p_actions[mod(action,2)+1]

                # ForwardLearner
                update!(learnerMB, state, action, next_state, reward)
                SPE_record[action_counter] = copy(learnerMB.spe[1])

                # Track V values (used as fmri regressor)
                V_state[j] = learnerMF.V[state]
                if done
                    j+=1
                    if j < (nr_samples + data_stats.nr_rewards)
                        V_state[j] = learnerMF.V[next_state] # Track V value of landing state in case it is the Goal
                    end
                end
                j+=1

                # Actor-Critic
                update!(learnerMF, state, action, next_state, reward, p_actions[action], modulatingfactor)
                TD_record[action_counter] = learnerMF.td_error[1]
            end # end t
        end # end i_episode
    end # end i_puzzle
    #---- Without bias
    trial_by_trial_p_sa = (action_selection_probs)
    if do_debug_print
        println(learnerHybrid.QSA)
        println("action_selection_probs: $action_selection_probs")
        println("LL_max: $LL")
        println("trial_by_trial_p_sa: $trial_by_trial_p_sa")
    end
    algo = "hybrid_actorcritic_FWD"

    param_v = [alpha, eta_critic, gamma, lambda_actor, lambda_critic, eta_fwd, w_MB, w_MF]
    param_str = "alpha=$(round(alpha,digits=4)),
                eta_critic=$(round(eta_critic,digits=4)),
                gamma=$(round(gamma,digits=4)),
                lambda_actor=$(round(lambda_actor,digits=4)),
                lambda_critic=$(round(lambda_critic,digits=4)),
                eta_fwd=$(round(eta_fwd,digits=4)),
                w_MB=$(round(w_MB,digits=4)),
                w_MF=$(round(w_MF,digits=4))"

    results_dict = Dict{String,Any}("finalQSA_Hybrid_record" => learnerHybrid.QSA,
                                    "final_PolicyParameters_MF" => learnerMF.policyparameters,
                                    "finalQSA_MB_record" => learnerMB.QSA,
                                    "final_Tsas" => learnerMB.TSAS,
                                    "trialByTrial_Vs" => V_state,
                                    "PolicyParametersDiffChosenVsUnchosen" => PolicyParametersDiffChosenVsUnchosen_record,
                                    "PolicyDiffChosenVsUnchosen" => PolicyDiffChosenVsUnchosen_record,
                                )
    signals = Step_by_step_signals(data_stats.all_subject_ids[1],
    TD_record, #RPE
    SPE_record, #SPE
    trial_by_trial_p_sa,
    LL,
    data_stats,
    algo,
    param_v,
    param_str,
    99, #biased_action_idx, dummy value
    results_dict)
    return LL, signals
end


text:
import DeltaRCWA

struct FunctionSheet{A <: Function, B <: Function} <: DeltaRCWA.Sheet
    η::A
    μ::B
end

DeltaRCWA.Zₑˣˣ(sheet::FunctionSheet, x) = 0.5sheet.μ(x[1])
DeltaRCWA.Yₘʸʸ(sheet::FunctionSheet, x) = 2sheet.η(x[1])

"""
    compute_DeltaRCWA_method(η::Function, μ::Function, L::Float64, N::Int, k::Float64, i::Int)

Compute the fields using DeltaRCWA (equivalent to Luke's RCWA method).
The output fields are on a grid compatible with the BIE method.

Arguments:
η:: A periodic function returning the metasurface magnetic conductance
μ:: A periodic function returning the metasurface electric resistance
L:: the length of the unit cell (some multiple of the period of η, μ)
N:: the number of discretization points along the unit cell
k:: the wavenumber/frequency of the incident planewave (assumed TM polarization)
i:: the index of the incident planewave in the Fourier basis

Returns:
x:: the positions of the grid points
kz₁:: the values of the wavevector in the z component for y<0
kz₂:: the values of the wavevector in the z component for y>0
Ĩ₁:: an array of the incident fields in the Fourier basis for y<0
Õ₁:: an array of the outgoing fields in the Fourier basis for y<0
Ĩ₂:: an array of the incident fields in the Fourier basis for y>0
Õ₂:: an array of the outgoing fields in the Fourier basis for y>0
"""
function compute_DeltaRCWA_method(η::Function, μ::Function, L::Float64, N::Int, k::Float64, i::Int)
    prob = DeltaRCWA.DeltaRCWAProblem(
        FunctionSheet(η, μ),
        ((N, L), ),
        k,
        [i == j ? 1.0 : 0.0 for j in 1:N],
        zeros(N),
    )
    sol = DeltaRCWA.solve(prob)
    x = sol.pw.x⃗[1]
    kz₁ = DeltaRCWA._get_kz(sol.pw, sol.stack.media[1])
    kz₂ = DeltaRCWA._get_kz(sol.pw, sol.stack.media[end])
    (x=x, kz₁=kz₁, kz₂=kz₂, Ĩ₁=sol.I₁, Õ₁=sol.O₁, Ĩ₂=sol.I₂, Õ₂=sol.O₂)
end


text:
using Statistics

"""
    peak_snr, peak_amp = peakSNR(snr, amp)

Estimate the snr peak by threading a line of slope 1 through the (`amp`,`snr`)
data. Both `amp` and `snr` are expressed in dB.
"""
function peakSNR(snr, amp)
    # Delete garbage data
    i = abs.(snr) .!= Inf
    snr = snr[i]
    amp = amp[i]
    i = .!isnan.(snr)
    snr = snr[i]
    amp = amp[i]
    i = snr .>= 3
    snr = snr[i]
    amp = amp[i]

    n = length(amp)

    # Sort by amplitude
    i = sortperm(amp)
    amp = amp[i]
    snr = snr[i]

    i = 1
    m = 0
    while any(i .!= 0) && n > 3
        # Draw a 45-degree median line through the data
        tmp = sort(snr-amp)
        if rem(n,2) == 0
            m = mean(tmp[n÷2 .+ [1,0]])
        else
            m = tmp[(n+1)÷2]
        end
        # Discard data that is more than 6dB off
        i = abs.(amp - snr .+ m) .> 6
        snr = snr[.!i]
        amp = amp[.!i]
        n = length(amp)
    end

    peak_amp = maximum(amp)
    peak_snr = peak_amp + m
    return peak_snr, peak_amp
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
export rule

@rule GCV(:y, Marginalisation) (m_x::Any, q_z::Any, q_κ::Any, q_ω::Any) = begin

    ksi = mean(q_κ) ^ 2 * var(q_z) + mean(q_z) ^ 2 * var(q_κ) + var(q_z) * cov(q_κ)
    A = exp(-mean(q_ω) + var(q_ω) / 2)
    B = exp(-mean(q_κ) * mean(q_z) + ksi / 2)

    return NormalMeanVariance(mean(m_x), var(m_x) + inv(A * B))
end

@rule GCV(:y, Marginalisation) (q_x::Any, q_z::Any, q_κ::Any, q_ω::Any) = begin

    ksi = mean(q_κ) ^ 2 * var(q_z) + mean(q_z) ^ 2 * var(q_κ) + var(q_z) * cov(q_κ)
    A = exp(-mean(q_ω) + var(q_ω) / 2)
    B = exp(-mean(q_κ) * mean(q_z) + ksi / 2)

    return NormalMeanVariance(mean(q_x), inv(A * B))
end


text:
@testset "TUDatasets" begin

    graphs = loadgraphs(TUDatasets.AIDSDataset())

    @test ng(graphs) == 2000
end


text:
module EquivariantNeuralOperators

include("./nn.jl")
export EquivConv,EquivAttn,EquivNL,LinearOperator

end # module

max_stars_repo_name: paulxshen/eno

max_issues_repo_name: paulxshen/eno

max_forks_repo_name: paulxshen/eno


text:
I = π * besseli(0.0, 1.0)
for endpt in [neither, left, right, both]
    x, w = chebyshev(10, 1, endpt)
    Q = sum(w .* exp.(x))
    @test I ≈ Q
end

I = π * besseli(1.0, 1.0)
for endpt in [neither, left, right, both]
    x, w = chebyshev(10, 2, endpt)
    Q = sum(w .* exp.(x))
    @test I ≈ Q
end


text:
@symbol_func function surface_area(cur_reactor::AbstractReactor)
  cur_left_part = 4.0

  cur_left_part *= cur_reactor.pi ^ 2

  cur_left_part *= a(cur_reactor)

  cur_left_part *= cur_reactor.R_0

  cur_right_part = ( cur_reactor.kappa_95 ^ 2 - 1 )

  cur_right_part *= ( 2 / cur_reactor.pi )

  cur_right_part += 1

  cur_right_part /= cur_reactor.kappa_95

  cur_area = cur_left_part

  cur_area *= cur_right_part

  cur_area
end


text:
function rho_kernel(
    n, k, Kf, Kf_mean, pq_end, eps_rp, eps_rp_min, rt_inc, rt_dec, eta,
    rho_max, rho_min_pq, rho_min_w,
    u_curr, u_prev, v_curr, v_prev, l_curr, l_prev, rho,
    tau, rp, rp_old, rp_k0
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if I <= n
        @inbounds begin
            delta_u = u_curr[I] - u_prev[I]
            delta_v = v_curr[I] - v_prev[I]
            delta_l = l_curr[I] - l_prev[I]
            alpha = abs(delta_l / delta_u)
            beta = abs(delta_l / delta_v)
            rho_v = rho[I]
            rp_v = rp[I]
            rp_old_v = rp_old[I]
            rp_k0_v = rp_k0[I]

            if abs(delta_l) <= eps_rp_min
                tau[k,I] = tau[k-1,I]
            elseif abs(delta_u) <= eps_rp_min && abs(delta_v) > eps_rp_min
                tau[k,I] = beta
            elseif abs(delta_u) > eps_rp_min && abs(delta_v) <= eps_rp_min
                tau[k,I] = alpha
            elseif abs(delta_u) <= eps_rp_min && abs(delta_v) <= eps_rp_min
                tau[k,I] = tau[k-1,I]
            else
                tau[k,I] = sqrt(alpha*beta)
            end

            if (k % Kf) == 0
                mean_tau = 0
                @inbounds for j=1:Kf_mean
                    mean_tau += tau[k-j+1,I]
                end
                mean_tau /= Kf_mean
                if mean_tau >= rt_inc*rho_v
                    if abs(rp_v) >= eps_rp && abs(rp_old_v) >= eps_rp
                        if abs(rp_v) > eta*abs(rp_k0_v) || abs(rp_old_v) > eta*abs(rp_k0_v)
                            rho_v *= rt_inc
                        end
                    end
                    #=
                elseif mean_tau > rt_inc2*rho_v
                    if abs(rp_v) >= eps_rp && abs(rp_old_v) >= eps_rp
                        if abs(rp_v) > eta*abs(rp_k0_v) || abs(rp_old_v) > eta*abs(rp_k0_v)
                            rho_v = rt_inc2*mean_tau
                        end
                    end
                    =#
                elseif mean_tau > rho_v
                    if abs(rp_v) >= eps_rp && abs(rp_old_v) >= eps_rp
                        if abs(rp_v) > eta*abs(rp_k0_v) || abs(rp_old_v) > eta*abs(rp_k0_v)
                            rho_v = mean_tau
                        end
                    end
                elseif mean_tau <= rho_v/rt_dec
                    rho_v /= rt_dec
                    #=
                elseif mean_tau < rho_v/rt_dec2
                    rho_v = mean_tau/rt_dec2
                    =#
                elseif mean_tau < rho_v
                    rho_v = mean_tau
                end
            end

            rho_v = min(rho_max, rho_v)
            if I <= pq_end
                rho_v = max(rho_min_pq, rho_v)
            else
                rho_v = max(rho_min_w, rho_v)
            end
            rho[I] = rho_v
        end
    end

    return
end

function rho_kernel_cpu(
    n, k, Kf, Kf_mean, pq_end, eps_rp, eps_rp_min, rt_inc, rt_dec, eta,
    rho_max, rho_min_pq, rho_min_w,
    u_curr, u_prev, v_curr, v_prev, l_curr, l_prev, rho,
    tau, rp, rp_old, rp_k0
)
    Threads.@threads for I=1:n
        @inbounds begin
            delta_u = u_curr[I] - u_prev[I]
            delta_v = v_curr[I] - v_prev[I]
            delta_l = l_curr[I] - l_prev[I]
            alpha = abs(delta_l / delta_u)
            beta = abs(delta_l / delta_v)
            rho_v = rho[I]
            rp_v = rp[I]
            rp_old_v = rp_old[I]
            rp_k0_v = rp_k0[I]

            if abs(delta_l) <= eps_rp_min
                tau[k,I] = tau[k-1,I]
            elseif abs(delta_u) <= eps_rp_min && abs(delta_v) > eps_rp_min
                tau[k,I] = beta
            elseif abs(delta_u) > eps_rp_min && abs(delta_v) <= eps_rp_min
                tau[k,I] = alpha
            elseif abs(delta_u) <= eps_rp_min && abs(delta_v) <= eps_rp_min
                tau[k,I] = tau[k-1,I]
            else
                tau[k,I] = sqrt(alpha*beta)
            end

            if (k % Kf) == 0
                mean_tau = mean(tau[k - Kf_mean+1:k,I])
                if mean_tau >= rt_inc*rho_v
                    if abs(rp_v) >= eps_rp && abs(rp_old_v) >= eps_rp
                        if abs(rp_v) > eta*abs(rp_k0_v) || abs(rp_old_v) > eta*abs(rp_k0_v)
                            rho_v *= rt_inc
                        end
                    end
                    #=
                elseif mean_tau > rt_inc2*rho_v
                    if abs(rp_v) >= eps_rp && abs(rp_old_v) >= eps_rp
                        if abs(rp_v) > eta*abs(rp_k0_v) || abs(rp_old_v) > eta*abs(rp_k0_v)
                            rho_v = rt_inc2*mean_tau
                        end
                    end
                    =#
                elseif mean_tau > rho_v
                    if abs(rp_v) >= eps_rp && abs(rp_old_v) >= eps_rp
                        if abs(rp_v) > eta*abs(rp_k0_v) || abs(rp_old_v) > eta*abs(rp_k0_v)
                            rho_v = mean_tau
                        end
                    end
                elseif mean_tau <= rho_v/rt_dec
                    rho_v /= rt_dec
                    #=
                elseif mean_tau < rho_v/rt_dec2
                    rho_v = mean_tau/rt_dec2
                    =#
                elseif mean_tau < rho_v
                    rho_v = mean_tau
                end
            end

            rho_v = min(rho_max, rho_v)
            if I <= pq_end
                rho_v = max(rho_min_pq, rho_v)
            else
                rho_v = max(rho_min_w, rho_v)
            end
            rho[I] = rho_v
        end
    end

    return
end


text:
using AccessMacro
using Base.Test

@access struct A
    x :: Int
end

@access mutable struct B{T}
    x
    y :: T
end

@access type C{S, T}
    x :: S
    y :: T
end

@access immutable D
    x
end

@test x(A(3)) == 3

@test x(B(1, "2")) == 1 && y(B(1, "2")) == "2"

@test x(C(1.0, [2])) == 1.0 && y(B(1.0, [2])) == [2]

@test x(D(43)) == 43


text:
using StateSpaceModels, BenchmarkTools, CSV
# Local level series
locallevel_series = CSV.read("./benchmark/locallevel.csv"; header = false)
verbose = 0

locallevel = local_level(Matrix{Float64}(locallevel_series))
locallevel.H[1] = 1
locallevel.Q[1] = 1
@benchmark StateSpaceModels.kalman_filter($locallevel)
# 24th October 2019
# BenchmarkTools.Trial:
#   memory estimate:  12.17 KiB
#   allocs estimate:  93
#   --------------
#   minimum time:     9.258 μs (0.00% GC)
#   median time:      9.832 μs (0.00% GC)
#   mean time:        12.296 μs (11.59% GC)
#   maximum time:     4.382 ms (99.48% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

locallevel = local_level(Matrix{Float64}(locallevel_series))
locallevel.H[1] = 1
locallevel.Q[1] = 1
@benchmark StateSpaceModels.univariate_kalman_filter($locallevel)
# 24th October 2019
# BenchmarkTools.Trial:
#   memory estimate:  8.33 KiB
#   allocs estimate:  36
#   --------------
#   minimum time:     4.977 μs (0.00% GC)
#   median time:      5.458 μs (0.00% GC)
#   mean time:        6.477 μs (10.23% GC)
#   maximum time:     406.721 μs (97.11% GC)
#   --------------
#   samples:          10000
#   evals/sample:     6

# Benchmark structural model
AP = CSV.read("./examples/AirPassengers.csv")
logAP = log.(Vector{Float64}(AP[:Passengers]))

structural_model = structural(logAP, 12)
structural_model.H[1] = 1e-4
structural_model.Q[[1, 5, 9]] .= 1e-4

@benchmark StateSpaceModels.kalman_filter($structural_model)
# 24th October 2019
# BenchmarkTools.Trial:
#   memory estimate:  740.64 KiB
#   allocs estimate:  1455
#   --------------
#   minimum time:     336.863 μs (0.00% GC)
#   median time:      352.551 μs (0.00% GC)
#   mean time:        473.981 μs (13.59% GC)
#   maximum time:     5.512 ms (77.69% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

structural_model = structural(logAP, 12)
structural_model.H[1] = 1e-4
structural_model.Q[[1, 5, 9]] .= 1e-4

@benchmark StateSpaceModels.univariate_kalman_filter($structural_model)
# 24th October 2019
# BenchmarkTools.Trial:
#   memory estimate:  658.19 KiB
#   allocs estimate:  446
#   --------------
#   minimum time:     246.795 μs (0.00% GC)
#   median time:      256.044 μs (0.00% GC)
#   mean time:        304.342 μs (12.20% GC)
#   maximum time:     3.262 ms (79.72% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1


text:
import Base: ==

# https://github.com/JuliaLang/julia/issues/11618
signed_prod(lst) = length(lst) == 1 ? lst[1] : Base.prod(lst)
signed_sum(lst) = length(lst) == 1 ? lst[1] : Base.sum(lst)

abstract Node
==(n1::Node, n2::Node) = isequal(n1.val, n2.val)
calc(n::Float64) = n
type Inv<:Node val end
calc(i::Inv) = 1.0/calc(i.val)
type Prd<:Node val end
calc(p::Prd) = signed_prod(map(calc, p.val))
type Neg<:Node val end
calc(n::Neg) = -calc(n.val)
type Sum<:Node val end
calc(s::Sum) = signed_sum(map(calc, s.val))

@with_names begin

    spc = Drop(Star(Space()))

    @with_pre spc begin

        sum = Delayed()
        val = E"(" + spc + sum + spc + E")" | PFloat64()

        neg = Delayed()             # allow multiple negations (eg ---3)
        neg.matcher = Nullable{Matcher}(val | (E"-" + neg > Neg))

        mul = E"*" + neg
        div = E"/" + neg > Inv
        prd = neg + (mul | div)[0:end] |> Prd

        add = E"+" + prd
        sub = E"-" + prd > Neg
        sum.matcher = Nullable{Matcher}(prd + (add | sub)[0:end] |> Sum)

        all = sum + spc + Eos()

    end
end

parse_one(" 1 + 2 * 3 / 4 ", Trace(all); debug=true)

for (src, val) in [
                   (" 1 ", 1),
                   (" - 1 ", -1),
                   (" 1 + 1 ", 2),
                   (" 1 - 1 ", 0),
                   (" - 1 - 1 ", -2)
                   ]
#    @test_approx_eq calc(parse_dbg(src, Trace(all))[1]) val
    @test_approx_eq calc(parse_one(src, Trace(all))[1]) val
    println("$src = $val")
end

println("fix ok")


text:
module BenchAvgFilter1D

import Random
using BenchmarkTools
using TapirBenchmarks

Random.seed!(1234)

const SUITE = BenchmarkGroup()

for n in [2^25, 2^27]
    SUITE[:n=>n] = s0 = BenchmarkGroup()
    ys0, xs = avgfilter1d_setup(n)
    s0[:impl=>:seq] = @benchmarkable avgfilter1d_seq!($(copy(ys0)), $xs)
    s0[:impl=>:threads] = @benchmarkable avgfilter1d_threads!($(copy(ys0)), $xs)
    s0[:impl=>:tapir_dac] = @benchmarkable avgfilter1d_tapir_dac!($(copy(ys0)), $xs)
    # s0[:impl=>:tapir_seq] = @benchmarkable avgfilter1d_tapir_seq!($(copy(ys0)), $xs)
end

end
BenchAvgFilter1D.SUITE


text:
"""
The sum of the squares of the first ten natural numbers is,
1^2 +2^2 +... +10^2 = 385
The square of the sum of the first ten natural numbers is,
(1+2+...+10)^2 =55^2 =3025
Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025-385 =2640
Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.
"""
function sumsqdiff(n)
sum(1:n)^2 - sum((1:n).^2)
end

sumsqdiff(100)


text:
using Yao, QuAlgorithmZoo
using Yao.ConstGate
using LinearAlgebra, Test, Random

@testset "BP diff" begin
    reg = rand_state(4)
    block = put(4, 2=>rot(X, 0.3))
    df = BPDiff(block)
    @test df.grad == 0
    @test nqubits(df) == 4

    df2 = BPDiff(rot(CNOT, 0.3))
    @test df2.grad == 0
    @test nqubits(df2) == 2
end

@testset "Qi diff" begin
    reg = rand_state(4)
    df2 = QDiff(rot(CNOT, 0.3))
    @test df2.grad == 0
    @test nqubits(df2) == 2

    @test df2' isa QDiff
    @test mat(df2) == mat(df2')'
end

"""
    loss_expect!(circuit::AbstractBlock, op::AbstractBlock) -> Function

Return function "loss!(ψ, θ) -> Vector"
"""
function loss_expect!(circuit::AbstractBlock, op::AbstractBlock)
    N = nqubits(circuit)
    function loss!(ψ::AbstractRegister, θ::Vector)
        params = parameters(circuit)
        dispatch!(circuit, θ)
        ψ |> circuit
        popdispatch!(circuit, params)
        expect(op, ψ)
    end
end

"""
    loss_Z1!(circuit::AbstractBlock; ibit::Int=1) -> Function

Return the loss function f = <Zi> (means measuring the ibit-th bit in computation basis).
"""
loss_Z1!(circuit::AbstractBlock; ibit::Int=1) = loss_expect!(circuit, put(nqubits(circuit), ibit=>Z))

cnot_entangler(n::Int, pairs) = chain(n, control(n, [ctrl], target=>X) for (ctrl, target) in pairs)

function rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false)
    rt = chain(nbit, [put(nbit, ibit=>Rz(0.0)), put(nbit, ibit=>Rx(0.0)), put(nbit, ibit=>Rz(0.0))])
    noleading && popfirst!(rt)
    notrailing && pop!(rt)
    rt
end

rset(nbit::Int, noleading::Bool=false, notrailing::Bool=false) = chain(nbit, [rotor(nbit, j, noleading, notrailing) for j=1:nbit])

function ibm_diff_circuit(nbit, nlayer, pairs)
    circuit = chain(nbit)

    ent = cnot_entangler(nbit, pairs)
    for i = 1:(nlayer + 1)
        i!=1 && push!(circuit, ent)
        push!(circuit, rset(nbit, i==1, i==nlayer+1))
    end
    circuit
end

@testset "BP diff" begin
    c = put(4, 3=>Rx(0.5)) |> autodiff(:BP)
    cad = c'
    @test mat(cad) == mat(c)'

    circuit = chain(4, repeat(4, H, 1:4), put(4, 3=>Rz(0.5)) |> autodiff(:BP), control(4, 2, 1=>shift(0.4)) |> autodiff(:BP), control(2, 1=>X), put(4, 4=>Ry(0.2)) |> autodiff(:BP))
    op = put(4, 3=>Y)
    θ = [0.1, 0.2, 0.3]
    dispatch!(circuit, θ)
    loss! = loss_expect!(circuit, op)
    ψ0 = rand_state(4)
    ψ = copy(ψ0) |> circuit

    # get gradient
    δ = ψ |> op
    backward!(δ, circuit)
    g1 = gradient(circuit)

    g2 = zero(θ)
    η = 0.01
    for i in 1:length(θ)
        θ1 = copy(θ)
        θ2 = copy(θ)
        θ1[i] -= 0.5η
        θ2[i] += 0.5η
        g2[i] = (loss!(copy(ψ0), θ2) - loss!(copy(ψ0), θ1))/η |> real
    end
    g3 = opdiff.(() -> copy(ψ0) |> circuit, collect_blocks(BPDiff, circuit), Ref(op))
    @test isapprox.(g1, g2, atol=1e-5) |> all
    @test isapprox.(g2, g3, atol=1e-5) |> all
end

@testset "constructor" begin
    @test generator(put(4, 1=>Rx(0.1))) == put(4, 1=>X)
    @test generator(Rx(0.1)) == X
    circuit = chain(put(4, 1=>Rx(0.1)), control(4, 2, 1=>Ry(0.3)))
    c2 = circuit |> autodiff(:BP)
    @test c2[1] isa BPDiff
    @test !(c2[2] isa BPDiff)
end

@testset "numdiff & opdiff" begin
    @test collect_blocks(XGate, chain([X, Y, Z])) == [X]

    c = chain(put(4, 1=>Rx(0.5))) |> autodiff(:QC)
    nd = numdiff(c[1].content) do
        expect(put(4, 1=>Z), zero_state(4) |> c) |> real # return loss please
    end

    ed = opdiff(c[1].content, put(4, 1=>Z)) do
        zero_state(4) |> c  # a function get output
    end
    @test isapprox(nd, ed, atol=1e-4)

    reg = rand_state(4)
    c = chain(put(4, 1=>Rx(0.5)), control(4, 1, 2=>Ry(0.5)), control(4, 1, 2=>shift(0.3)),  kron(4, 2=>Rz(0.3), 3=>Rx(0.7))) |> autodiff(:QC)
    dbs = collect_blocks(QDiff, c)
    loss1z() = expect(kron(4, 1=>Z, 2=>X), copy(reg) |> c) |> real  # return loss please
    nd = numdiff.(loss1z, dbs)
    ed = opdiff.(()->copy(reg) |> c, dbs, Ref(kron(4, 1=>Z, 2=>X)))
    gd = gradient(c)
    @test gradient(c, :QC) == gd
    @test gradient(c, :BP) == []
    @test isapprox(nd, ed, atol=1e-4)
    @test ed == gd
end

@testset "stat" begin
    nbit = 3
    f(x::Number, y::Number) = Float64(abs(x-y) < 1.5)
    x = 0:1<<nbit-1
    h = f.(x', x)
    V = StatFunctional(h)
    VF = StatFunctional{2}(f)
    prs = [1=>2, 2=>3, 3=>1]
    c = ibm_diff_circuit(nbit, 2, prs) |> autodiff(:QC)
    dispatch!(c, :random)
    dbs = collect_blocks(AbstractDiff,c)

    p0 = zero_state(nbit) |> c |> probs
    sample0 = measure(zero_state(nbit) |> c; nshots=5000)
    loss0 = expect(V, p0 |> as_weights)
    gradsn = numdiff.(()->expect(V, zero_state(nbit) |> c |> probs |> as_weights), dbs)
    gradse = statdiff.(()->zero_state(nbit) |> c |> probs |> as_weights, dbs, Ref(V), initial=p0 |> as_weights)
    gradsf = statdiff.(()->measure(zero_state(nbit) |> c; nshots=5000), dbs, Ref(VF), initial=sample0)
    @test all(isapprox.(gradse, gradsn, atol=1e-4))
    @test norm(gradsf-gradse)/norm(gradsf) <= 0.2

    # 1D
    h = randn(1<<nbit)
    V = StatFunctional(h)
    c = ibm_diff_circuit(nbit, 2, prs) |> autodiff(:QC)
    dispatch!(c, :random)
    dbs = collect_blocks(AbstractDiff, c)

    p0 = zero_state(nbit) |> c |> probs |> as_weights
    loss0 = expect(V, p0 |> as_weights)
    gradsn = numdiff.(()->expect(V, zero_state(nbit) |> c |> probs |> as_weights), dbs)
    gradse = statdiff.(()->zero_state(nbit) |> c |> probs |> as_weights, dbs, Ref(V))
    @test all(isapprox.(gradse, gradsn, atol=1e-4))
end


text:
@safetestset "SVG API support in Renderer.Html" begin

  @safetestset "" begin
    using Genie
    using Genie.Renderer.Html

    @test_throws UndefVarError svg()
  end;

  @safetestset "Update config when custom startup args" begin
    using Genie
    using Genie.Renderer.Html

    Html.register_svg_slements()
    @test svg() == "<svg></svg>"

    @test_throws UndefVarError clippath()
    @test clipPath() == "<clipPath></clipPath>"
  end;

end;


text:
const PARALLEL = false
macro ifparallel(e)
    if PARALLEL
        return e
    end
end
@ifparallel if length(workers()) == 1
    addprocs(4)
end
pyplot()
default(grid=false)
# @everywhere using Revise
@ifparallel using ParallelDataTransfer
@everywhere using Parameters, JacProp, OrdinaryDiffEq, LTVModels, LTVModelsBase
@everywhere using Flux: params, jacobian
@everywhere using Flux.Optimise: Param, optimiser, expdecay
@everywhere begin

    @with_kw struct LinearSys <: AbstractSystem
        A
        B
        N = 1000
        nx = size(A,1)
        nu = size(B,2)
        h = 0.02
        σ0 = 0
        sind = 1:nx
        uind = nx+1:(nx+nu)
        s1ind = (nx+nu+1):(nx+nu+nx)
    end

    function LinearSys(seed; nx = 10, nu = nx, h=0.02, kwargs...)
        srand(seed)
        A = randn(nx,nx)
        A = A-A'        # skew-symmetric = pure imaginary eigenvalues
        A = A - h*I     # Make 'slightly' stable
        A = expm(h*A)   # discrete time
        B = h*randn(nx,nu)
        LinearSys(;A=A, B=B, nx=nx, nu=nu, h=h, kwargs...)
    end

    function generate_data(sys::LinearSys, seed, validation=false)
        Parameters.@unpack A,B,N, nx, nu, h, σ0 = sys

        srand(seed)

        u      = filt(ones(5),[5], 10randn(N+2,nu))'
        t      = h:h:N*h+h
        x0     = randn(nx)
        x      = zeros(nx,N+1)
        x[:,1] = x0
        for i = 1:N-1
            x[:,i+1] = A*x[:,i] + B*u[:,i]
        end

        validation || (x .+= σ0 * randn(size(x)))
        u = u[:,1:N]
        @assert all(isfinite, u)
        x,u
    end

    function true_jacobian(sys::LinearSys, x, u)
        [sys.A sys.B]
    end

    lastprint = 0
    function callbacker(epoch, loss,d,trace,model)
        function ()
            global lastprint
            Flux.reset!(model.m)
            l = sum(d->Flux.data(loss(d...)),d)
            Flux.reset!(model.m)
            increment!(trace,epoch,l)
            if epoch % 20 == 0 && lastprint != epoch
                @show epoch
                println(@sprintf("Loss: %.4f", l))
                plot(trace, reuse=true, yscale=:log10)
                gui()
                lastprint = epoch
            end
        end
    end

    num_params = 20
    wdecay     = 0
    stepsize   = 0.005
    const sys  = LinearSys(1, N=200, h=0.02, σ0 = 1)
    true_jacobian(x,u) = true_jacobian(sys,x,u)
    nu         = sys.nu
    nx         = sys.nx

end

numtrajs = 3
vt = Trajectory(generate_data(sys, 100)...) #' Generate validation data
trajs = [Trajectory(generate_data(sys, i)...) for i = 1:numtrajs] #' Generate training data
@ifparallel sendto(collect(2:nprocs()), trajs=trajs, vt=vt) # Send data to workers

# f1 = @spawnat 4 begin
srand(1)
models     = [RecurrentSystem(nx,nu,num_params, a) for a in default_activations]
opts       = ADAM.(params.(models), stepsize, decay=0.0005)
trainers  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker)
for i = 1:numtrajs
    trainers(trajs[i], epochs=0, jacprop=0)
end
trainers
# endtrainerswd
trainers(epochs=2000, jacprop=0)
predsimplot(trainers.models, vt, layout=10)
traceplot(trainers, reuse=false)

# f2 = @spawnat 2 begin
srand(1)
models     = [RecurrentDiffSystem(nx,nu,num_params, a) for a in default_activations]
opts       = ADAM.(params.(models), stepsize, decay=0.0005)
trainerds  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker)
for i = 1:numtrajs
    trainerds(trajs[i], epochs=0, jacprop=0)
end
trainerds
# end
trainerds(epochs=2000, jacprop=0)
predsimplot(trainerds.models, vt, layout=10)
traceplot(trainerds, reuse=false)

wdecay = 0.01
# f3 = @spawnat 3 begin
srand(1)
models     = [RecurrentSystem(nx,nu,num_params, a) for a in default_activations]
opts       = [[ADAM(params(models[i]), stepsize, decay=0.0005); [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
trainerswd  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker)
for i = 1:numtrajs
    trainerswd(trajs[i], epochs=0, jacprop=0)
end
trainerswd
# end
trainerswd(epochs=2000, jacprop=0)
predsimplot(trainerswd.models, vt, layout=10)
traceplot(trainerswd, reuse=false)

# f4 = @spawnat 1 begin
srand(1)
models     = [RecurrentDiffSystem(nx,nu,num_params, a) for a in default_activations[3:3]]
opts       = [[ADAM(params(models[i]), stepsize, decay=0.0005); [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
trainerdswd  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker)
for i = 1:numtrajs
    trainerdswd(trajs[i], epochs=0, jacprop=0)
    # traceplot(trainerdswd)
end
trainerdswd
# end
trainerdswd(epochs=20000, jacprop=0)
predsimplot(trainerdswd.models, vt, layout=10)
traceplot(trainerdswd, reuse=false)

# trainers,trainerds = fetch(f1), fetch(f2)
# trainerswd,trainerdswd = fetch(f3), fetch(f4)

traceplot(trainers, layout=4)
traceplot!(trainerds, subplot=2)
traceplot!(trainerswd, subplot=3)
traceplot!(trainerdswd, subplot=4)

eigvalplot(trainers.models, vt, true_jacobian, title="f, wdecay=0", layout=4, subplot=1, markeralpha=0.05, ds=1, markerstrokealpha=0)
eigvalplot!(trainerds.models, vt, true_jacobian,  title="g, wdecay=0", subplot=2, markeralpha=0.05, ds=1, markerstrokealpha=0)
eigvalplot!(trainerswd.models, vt, true_jacobian,  title="f, wdecay=$wdecay", subplot=3, markeralpha=0.05, ds=1, markerstrokealpha=0)
eigvalplot!(trainerdswd.models, vt, true_jacobian,  title="g, wdecay=$wdecay", subplot=4, markeralpha=0.05, ds=1, markerstrokealpha=0)
plot!(link=:both, xlims=(0,1.1), ylims=(-0.5, 0.5))
gui()

ui = display_modeltrainer(trainerdswd, size=(800,600))

for acti in 1:4
    plot(layout=4, link=:both)
    actstr = string(JacProp.default_activations[acti])
    actstr = actstr[1:2] == "NN" ? actstr[7:end] : actstr
    for tri in 1:4
        tr = [trainers, trainerds, trainerswd, trainerdswd][tri]
        eigvalplot!([tr.models[acti]], vt, true_jacobian, title="$actstr, $(tri%2==0 ? "g" : "f"), wdecay=$(tri>2 ? wdecay : 0)", subplot=tri, markeralpha=0.05, ds=3, markerstrokealpha=0)
    end
    gui()
end

JLD.@save "workspace" trainers trainerds trainerswd trainerdswd trajs


text:
@testset "bayesian_gev.jl" begin
    n = 100

    x₁ = Variable("x₁", randn(n))
    x₂ = Variable("x₂", randn(n) / 3)
    x₃ = Variable("x₃", randn(n) / 10)

    μ = 1.0 .+ x₁.value
    ϕ = -.05 .+ x₂.value
    ξ = x₃.value

    σ = exp.(ϕ)
    θ = [1.0; 1.0; -0.05; 1.0; 0.0; 1.0]

    pd = GeneralizedExtremeValue.(μ, σ, ξ)
    y = rand.(pd)

    @testset "gevfitbayes(y; locationcov, logscalecov, shapecov, niter, warmup)" begin
        # model building with non-stationary location, logscale and shape
        fm = Extremes.gevfitbayes(y,
            locationcov = [x₁],
            logscalecov = [x₂],
            shapecov = [x₃],
            niter=500, warmup=5)

        # data is y
        @test fm.model.data.value ≈ y

        # location is x₁
        @test length(fm.model.location.covariate) == 1
        @test fm.model.location.covariate[1].value ≈ x₁.value

        # logscale is x₂
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₂.value

        # shape is x₃
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₃.value

    end

    @testset "gevfitbayes(df, datacol; locationcovid, logscalecovid, shapecovid, niter, warmup)" begin
        # model building with non-stationary location, logscale and shape
        df = DataFrame(y = y, x1 = x₁.value, x2 = x₂.value, x3 = x₃.value)

        fm = Extremes.gevfitbayes(df, :y, locationcovid = [:x1], logscalecovid = [:x2], shapecovid = [:x3], niter=500, warmup=5)

        # data is y
        @test fm.model.data.value ≈ y

        # location is x₁
        @test length(fm.model.location.covariate) == 1
        @test fm.model.location.covariate[1].value ≈ x₁.value

        # logscale is x₂
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₂.value

        # shape is x₃
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₃.value

    end

    @testset "gevfitbayes(model; niter, warmup)" begin
        # non-stationary location, logscale and shape
        model = Extremes.BlockMaxima(Variable("y", y),
            locationcov = [x₁],
            logscalecov = [x₂],
            shapecov = [x₃])

        fm = Extremes.gevfitbayes(model, niter=500, warmup=5)

        # data is y
        @test fm.model.data.value ≈ y

        # location is x₁
        @test length(fm.model.location.covariate) == 1
        @test fm.model.location.covariate[1].value ≈ x₁.value

        # logscale is x₂
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₂.value

        # shape is x₃
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₃.value

    end

end