text:
using PrincipalMomentAnalysis
using Test
using LinearAlgebra
using Statistics

# change back to this simple oneliner if we can just ignore problems in StaticArrays that is imported by NearestNeighbors
# @test [] == detect_ambiguities(Base, Core, PrincipalMomentAnalysis)
@testset "Method ambiguities" begin
	ambiguities = detect_ambiguities(Base, Core, PrincipalMomentAnalysis)
	filter!(x->x[1].module==PrincipalMomentAnalysis || x[2].module==PrincipalMomentAnalysis, ambiguities)
	@test ambiguities == []
end

const simplices2kernelmatrix     = PrincipalMomentAnalysis.simplices2kernelmatrix
const simplices2kernelmatrixroot = PrincipalMomentAnalysis.simplices2kernelmatrixroot

include("test_simplices.jl")
include("test_pma.jl")
include("test_util.jl")


text:
#!/usr/bin/env julia

using Plots
using VectorizedRoutines

plotly()  # Select a plotting backend.

# Read the antenna data file (CSV format).
raw_data = readdlm("dipole_pattern.csv", ',', Float32, skipstart=1)

# Separate and sort the data.
freq = unique(raw_data[:,1])
theta = unique(raw_data[:,2])
phi = unique(raw_data[:,3])
pattern_db = reshape(raw_data[:,4], (length(phi), length(theta), length(freq)))

# If complex gain values are needed, uncomment the following line:
#complex_pattern = 10^(raw_data[:,4]/10)*exp(1j*raw_data[:,5]/180*π)

## Plot in 2D ##

# Pick a frequency
fidx = 1

# Control plot dynamic range
pltmax = 5.
pltmin = -35.

# Create the 2D colormap (heatmap)
xs = [string(i) for i = round(Int32, theta)]
ys = [string(i) for i = round(Int32, phi)]
heatmap(ys, xs, pattern_db[:,:,fidx]', clims = (pltmin,pltmax),
        title = "Dipole Realized Gain", xlabel = "phi (deg)",
        ylabel = "theta (deg)", fc=:haline, colorbar_title="dBi")

## Plot in 3D (project onto surface of sphere) ##
theta2d, phi2d = Matlab.meshgrid(theta, phi)
x2d = sind(theta2d).*cosd(phi2d)
y2d = sind(theta2d).*sind(phi2d)
z2d = cosd(theta2d)
color_mat = pattern_db[:,:,fidx]
color_mat[pattern_db[:,:,fidx] .< pltmin] = pltmin

surface(x2d, y2d, z2d, fill_z = color_mat, fc = :haline, grid=false,
        title = "Dipole Realized Gain", colorbar_title="dBi")

# Remove the axes
# Currently ticks=false doesn't entirely remove the axes.
plot!(foreground_color = :white, foreground_color_title = :black)

## Plot in 3D ##
r = pattern_db[:,:,fidx] - pltmin
r[r .< 0.0] = 0.0

x2d = r.*sind(theta2d).*cosd(phi2d)
y2d = r.*sind(theta2d).*sind(phi2d)
z2d = r.*cosd(theta2d)

surface(x2d, y2d, z2d, fill_z = color_mat, fc = :haline,
        title = "Dipole Realized Gain", colorbar_title="dBi",
        grid = false)

# Remove the axes
# Currently ticks=false doesn't entirely remove the axes.
plot!(foreground_color = :white, foreground_color_title = :black)

ext: jl


text:
using Distributed
using  Hwloc
if length(workers()) == 1
    addprocs(Hwloc.num_physical_cores())  # get ready for parallelism
end
@everywhere begin
    using HypothesisTests: pvalue, UnequalVarianceTTest, VarianceFTest, EqualVarianceTTest
    using MLJ
    using DataFrames: DataFrame
    using CategoricalArrays: CategoricalArray
    using Distributions: UnivariateDistribution
end
include("ICPlibrary.jl")
@everywhere MLJ.@load LinearRegressor          pkg = GLM
@everywhere MLJ.@load LinearBinaryClassifier   pkg = GLM

export getPValueLinear

@everywhere function getPValueLinear(
    X::DataFrame,   # X must have MLJ scientific types and hot encoded already
    y::Vector{Float64},
    ExpInd::DataFrame,
    numberOfTargets::Int64,
    numberOfEnvironments::Int64,
)
    println("getPValueLinear: ", names(X))
    if numberOfTargets > 2
        LinearRegressorPipe = @pipeline Standardizer LinearRegressor
        LinearModel = machine(LinearRegressorPipe, X, y)
        fit!(LinearModel)
        coefs = fitted_params(LinearModel).linear_regressor.coef
        ŷ = MLJ.predict(LinearModel, X)
        yhatResponse = [ŷ[i,1].μ for i in 1:length(ŷ)]
        residuals = y .- yhatResponse
        r = report(LinearModel)
        coefsStdError = r.linear_regressor.stderror[2:end]  #### is the fist the intercept??
    else
        LinearBinaryClassifierPipe = @pipeline Standardizer LinearBinaryClassifier
        yc =  categorical(y[:, 1])
        LogisticModel = machine(LinearBinaryClassifierPipe, X, yc)
        fit!(LogisticModel)
        fp = fitted_params(LogisticModel)
        k = collect(keys(fp.fitted_params_given_machine))[2]
        coefs = fp.fitted_params_given_machine[k].coef
        ŷ = MLJ.predict(LogisticModel, X)
        yhatResponse = [pdf(ŷ[i], 1.0) for i in 1:length(ŷ)]     # probability to 1 ala R
        residuals = y - yhatResponse
        r = report(LogisticModel)
        k = collect(keys(r.report_given_machine))[2]
        coefsStdError = r.report_given_machine[k].stderror[2:end]
    end

    pvalVector = Vector{Float64}(undef, numberOfEnvironments)
    if numberOfTargets > 2
        for e in unique(ExpInd[:, 1])
            pvalVector[e] = pvalDoubler(
                getResiduals(ExpInd[:, 1], e, residuals, "inEnvironment"),
                getResiduals(ExpInd[:, 1], e, residuals, "notinEnvironment"),
            )
        end
    else
        for e in unique(ExpInd[:, 1])
            pvalVector[e] = getpvalClassification(
                getResiduals(ExpInd[:, 1], e, residuals, "inEnvironment"),
                getResiduals(ExpInd[:, 1], e, residuals, "notinEnvironment"),
            )
        end
    end

    pval = minimum(pvalVector) * (numberOfEnvironments - 1)
    pval = min(1, pval)
    return ((pvalue = pval, coefs = coefs, coefsStdError = coefsStdError))
end

@everywhere function getpvalClassification(
    xSample::Vector{Float64},
    ySample::Vector{Float64};
    nsim::Int64 = 500,
)
    return (pvalue(UnequalVarianceTTest(xSample, ySample)))
end # getpvalClassification

@everywhere function pvalDoubler(
    xSample::Vector{Float64},
    ySample::Vector{Float64};
)
    ttest = pvalue(UnequalVarianceTTest(xSample, ySample))
    vartest = pvalue(VarianceFTest(xSample, ySample))
    pval = 2 * min(ttest, vartest)

    return (pval)
end


text:
using Test
using GeoParams

@testset "EnergyParameters.jl" begin

    # This tests the MaterialParameters structure
    CharUnits_GEO   =   GEO_units(viscosity=1e19, length=10km);

    # Heat capacity ---------

    # Constant heat capacity
    cp1      =  ConstantHeatCapacity()
    info     =  param_info(cp1)
    @test isbits(cp1)
    @test cp1.cp.val == 1050.0

    cp1     = nondimensionalize(cp1,CharUnits_GEO)
    @test cp1.cp.val ≈ 1.3368075000000002e22
    @test compute_heatcapacity(cp1) ≈ 1.3368075000000002e22  # compute
    @test cp1() ≈ 1.3368075000000002e22  # compute
    @test compute_heatcapacity(cp1, random_name=1) ≈ 1.3368075000000002e22  # compute

    # Temperature-dependent heat capacity
    # dimensional
    T        =   250.0:100:1250;
    cp2      =   T_HeatCapacity_Whittington()
    Cp       =   similar(T)
    @test isbits(cp2)
    args = (;T=T)
    compute_heatcapacity!(Cp, cp2, args)
    @test sum(Cp) ≈ 11667.035717418683

    # nondimensional
    cp2_nd   =   T_HeatCapacity_Whittington()
    cp2_nd   =   nondimensionalize(cp2_nd,CharUnits_GEO)
    T_nd     =   Float64.(T*K/CharUnits_GEO.Temperature)
    Cp_nd    =   similar(T)
    args = (;T=T_nd)
    compute_heatcapacity!(Cp_nd, cp2_nd, args)
    @test sum(Cp_nd) ≈ 1.4853886523631602e23

    # Dimensionalize again and double-check the results
    @test sum(abs.(ustrip.(Cp_nd*CharUnits_GEO.heatcapacity) - Cp)) < 1e-11

    # Test with arrays
    T_array     =  T*ones(10)'
    Cp_array    =  similar(T_array)
    compute_heatcapacity!(Cp_array, cp1, (;))
    @test Cp_array[1] ≈ 1.3368075000000002e22

    Cp_array    =  similar(T_array)
    compute_heatcapacity!(Cp_array, cp2, (;T=T_array))
    @test sum(Cp_array[:,1]) ≈ 11667.035717418683

    T_array     =  T*ones(10)'
    Cp_array    =  zeros(size(T_array))
    compute_heatcapacity!(Cp_array, cp2, (;T=T_array))
    @test sum(Cp_array[:,1]) ≈ 11667.035717418683

    # Check that it works if we give a phase array
    MatParam    =   Array{MaterialParams, 1}(undef, 2);
    MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=1,
                        HeatCapacity  = ConstantHeatCapacity());

    MatParam[2] =   SetMaterialParams(Name="Crust", Phase=2,
                        HeatCapacity  = T_HeatCapacity_Whittington());

    Mat_tup = Tuple(MatParam)

    Mat_tup1 =   (  SetMaterialParams(Name="Mantle", Phase=1,
                        HeatCapacity  = ConstantHeatCapacity()),
                    SetMaterialParams(Name="Crust", Phase=2,
                        HeatCapacity  = ConstantHeatCapacity(cp=1100J/kg/K))
                 );

    # test computing material properties
    n = 100;
    Phases              = ones(Int64,n,n,n);
    Phases[:,:,20:end] .= 2;

    Cp = zeros(size(Phases));
    T  =  ones(size(Phases))*1500;
    P  =  zeros(size(Phases));

    args=(;T=T);
    compute_heatcapacity!(Cp, Mat_tup, Phases, args)    # computation routine w/out P (not used in most heat capacity formulations)
    @test sum(Cp[1,1,:]) ≈ 121399.0486067196

    # check with array of constant properties (and no required input args)
    args1=(;);
    compute_heatcapacity!(Cp, Mat_tup1, Phases, args1)    # computation routine w/out P (not used in most heat capacity formulations)
    @test sum(Cp[1,1,:]) ≈ 109050.0

    num_alloc = @allocated compute_heatcapacity!(Cp, Mat_tup, Phases, args)
    @test sum(Cp[1,1,:]) ≈ 121399.0486067196
    @test num_alloc <= 32

    # test if we provide phase ratios
    PhaseRatio  = zeros(n,n,n,3);
    for i in CartesianIndices(Phases)
        iz = Phases[i]
        I = CartesianIndex(i,iz)
        PhaseRatio[I] = 1.0
    end
    compute_heatcapacity!(Cp, Mat_tup, PhaseRatio, args)
    num_alloc = @allocated compute_heatcapacity!(Cp, Mat_tup, PhaseRatio, args)
    @test sum(Cp[1,1,:]) ≈ 121399.0486067196
    @test num_alloc <= 32

    # -----------------------

    # Conductivity ----------

    # Constant

    # Constant conductivity
    cond      =   ConstantConductivity()
    @test isbits(cond)
    @test NumValue(cond.k) == 3.0
    @test cond.k.unit==u"W"/K/m

    cond = nondimensionalize(cond,CharUnits_GEO)
    @test  NumValue(cond.k) ≈ 3.8194500000000007

    @test compute_conductivity(cond, T=100.0) ≈ 3.8194500000000007 # compute

    # Temperature-dependent conductivity
    # dimensional
    T        =   Vector{Float64}(250:100:1250);
    cond2    =   T_Conductivity_Whittington()
    k        =   compute_conductivity(cond2, T)
    @test isbits(cond2)
    @test sum(k) ≈ 27.503366436682285

    # nondimensional
    cond2_nd =   T_Conductivity_Whittington()
    cond2_nd =   nondimensionalize(cond2_nd,CharUnits_GEO)
    T_nd     =   Float64.(ustrip.(T/CharUnits_GEO.Temperature))
    k_nd     =   compute_conductivity(cond2_nd, T_nd)
    @test sum(k_nd) ≈ 35.01591097886205

    # Dimensionalize again and double-check the results
    @test sum(abs.(ustrip.(k_nd*CharUnits_GEO.conductivity) - k)) < 1e-11

    # Temperature-dependent parameterised conductivity
    # dimensional
    T        =   Vector{Float64}(250:100:1250);
    cond2    =   T_Conductivity_Whittington_parameterised()
    k        =   compute_conductivity(cond2, T)
    @test isbits(cond2)
    @test sum(k) ≈ 27.553653387829254

    # nondimensional
    cond2_nd =   T_Conductivity_Whittington_parameterised()
    cond2_nd =   nondimensionalize(cond2_nd,CharUnits_GEO)
    T_nd     =   Float64.(ustrip.(T/CharUnits_GEO.Temperature))
    k_nd     =   compute_conductivity(cond2_nd, T_nd)
    @test sum(k_nd) ≈ 35.079933810714806

    # Dimensionalize again and double-check the results
    @test sum(abs.(ustrip.(k_nd*CharUnits_GEO.conductivity) - k)) < 1e-11

    # Check if we use arrays
    T_array     =  ustrip.(T)*ones(100)'
    k_array     =  copy(T_array)
    P_array     =  copy(T_array)

    compute_conductivity!(k_array,cond)
    @test k_array[1] ≈ 3.8194500000000007

    args=(T=T_array, P=P_array)
    compute_conductivity!(k_array,cond, args)
    @test k_array[1] ≈ 3.8194500000000007

    compute_conductivity!(k_array,cond2, T=T_array)
    @test sum(k_array) ≈ 2755.3653387829254

    compute_conductivity!(k_array,cond2, (;T=T_array))
    @test sum(k_array) ≈ 2755.3653387829254

    k_TP    =   Set_TP_Conductivity("LowerCrust")
    compute_conductivity!(k_array, k_TP, (P=P_array, T=T_array))
    @test sum(k_array) ≈ 2055.7129327367625

    # Check that it works if we give a phase array
    MatParam    =   Array{MaterialParams, 1}(undef, 3);
    MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=1,
                        Conductivity  = ConstantConductivity());

    MatParam[2] =   SetMaterialParams(Name="Crust", Phase=2,
                        Conductivity  = T_Conductivity_Whittington());

    MatParam[3] =   SetMaterialParams(Name="MantleLithosphere", Phase=3,
                        Conductivity = Set_TP_Conductivity("Mantle"));

    Mat_tup = Tuple(MatParam)

    # test computing material properties
    n = 100;
    Phases              = ones(Int64,n,n,n);
    Phases[:,:,20:end] .= 2
    Phases[:,:,60:end] .= 3

    PhaseRatio  = zeros(n,n,n,3);
    for i in CartesianIndices(Phases)
        iz = Phases[i]
        I = CartesianIndex(i,iz)
        PhaseRatio[I] = 1.0
    end

    k   =  zeros(size(Phases))
    T   =  ones(size(Phases))*1500
    P   =  zeros(size(Phases))
    args= (P=P, T=T)

    compute_conductivity!(k, Mat_tup, Phases, args)
    num_alloc = @allocated compute_conductivity!(k, Mat_tup, Phases, args)
    @test sum(k) ≈ 1.9216938849389635e6
    @test num_alloc <= 32

    compute_conductivity!(k, Mat_tup, PhaseRatio, args)
    num_alloc = @allocated compute_conductivity!(k, Mat_tup, PhaseRatio, args)
    @test num_alloc <= 32
    @test sum(k) ≈ 1.9216938849389635e6

    #

    ######

    # TP-dependent conductivity for different predefines cases
    T       =   Vector{Float64}(250:100:1250);
    P       = 1e6*ones(size(T))/ustrip(uconvert(Pa,1MPa))  # must be in MPa!
    List    = ["LowerCrust"   "Mantle"        "OceanicCrust"  "UpperCrust"]
    Sol_kT  = [20.55712932736763 28.700405819019323 20.55712932736763 19.940302462417037]
    for i=1:length(List)
        k_TP    =   Set_TP_Conductivity(List[i])
        k       =   compute_conductivity(k_TP, P, T)         # note that P must be in MPa
        @test sum(k) ≈ Sol_kT[i]

        k_TP_nd  =   deepcopy(k_TP)
        k_TP_nd  =   nondimensionalize(k_TP_nd,CharUnits_GEO)
        T_nd     =   Float64.(ustrip.(T/CharUnits_GEO.Temperature))
        P_nd     =   Float64.(ustrip(P/CharUnits_GEO.stress))
        k_nd     =   compute_conductivity(k_TP_nd, P_nd, T_nd)

        @test ustrip(sum(abs.(ustrip.(k_nd*CharUnits_GEO.conductivity) - k))) < 1e-11

    end

    T = [200. 300.; 400. 500.]
    k1        =   compute_conductivity(cond2,T)
    # -----------------------

    # Latent heat -----------
    a = ConstantLatentHeat()
    Q_L = compute_latent_heat(a)
    @test isbits(a)
    @test Q_L == 400

    a   = nondimensionalize(a,CharUnits_GEO)
    Q_L = compute_latent_heat(a)
    @test Q_L ≈ 4e21

    # Check that it works if we give a phase array (including with an empty field)
    Mat_tup    =  ( SetMaterialParams(Name="Mantle", Phase=1,
                    LatentHeat  = ConstantLatentHeat()),
                SetMaterialParams(Name="Crust", Phase=2,
                    LatentHeat  = ConstantLatentHeat(Q_L=153kJ/kg)),
                SetMaterialParams(Name="MantleLithosphere", Phase=3)
                )

    # test computing material properties
    n = 100;
    Phases              = ones(Int64,n,n,n);
    Phases[:,:,20:end] .= 2;
    Phases[:,:,60:end] .= 3;

    PhaseRatio  = zeros(n,n,n,3);
    for i in CartesianIndices(Phases)
        iz = Phases[i]
        I = CartesianIndex(i,iz)
        PhaseRatio[I] = 1.0
    end

    Hl   =  zeros(size(Phases));
    z    =  ones(size(Phases))*10e3;
    args = (;)

    compute_latent_heat!(Hl, Mat_tup, Phases, args)
    @test  minimum(Hl)  ≈ 0.0
    @test  maximum(Hl)  ≈ 400
    @test  Hl[50,50,50] ≈ 153.0

    compute_latent_heat!(Hl, Mat_tup, PhaseRatio, args)
    @test sum(Hl)≈ 1.372e8

    # -----------------------

    # Radioactive heat ------
    a = ConstantRadioactiveHeat()
    H_r = compute_radioactive_heat(a)
    @test isbits(a)
    @test H_r ≈ 1.0e-6

    a = nondimensionalize(a,CharUnits_GEO)
    H_r = compute_radioactive_heat(a)
    @test H_r ≈ 0.1

    # depth-dependent radioactive heating:
    a = ExpDepthDependentRadioactiveHeat()
    z = 10e3
    H_r = compute_radioactive_heat(a,(z=z))
    @test isbits(a)
    @test H_r ≈3.678794411714423e-7

    Nx,Nz = 101,101
    z   = ones(Nx,Nz)*10e3;
    Hr=zero(z)
    compute_radioactive_heat!(Hr,a,(z=z))

    @test  sum(H_r) ≈ 3.678794411714423e-7

    # Check that it works if we give a phase array (including with an empty field)
    Mat_tup    =  ( SetMaterialParams(Name="Mantle", Phase=1,
                    RadioactiveHeat  = ConstantRadioactiveHeat()),
                SetMaterialParams(Name="Crust", Phase=2,
                    RadioactiveHeat  = ExpDepthDependentRadioactiveHeat()),
                SetMaterialParams(Name="MantleLithosphere", Phase=3)
                )
    Mat_tup1    =  ( SetMaterialParams(Name="Mantle", Phase=1,
                RadioactiveHeat  = ConstantRadioactiveHeat()),
            SetMaterialParams(Name="Crust", Phase=2,
                RadioactiveHeat  = ConstantRadioactiveHeat()),
            SetMaterialParams(Name="MantleLithosphere", Phase=3)
            )

    # test computing material properties
    n = 100;
    Phases              = ones(Int64,n,n,n);
    Phases[:,:,20:end] .= 2;
    Phases[:,:,60:end] .= 3;

    PhaseRatio  = zeros(n,n,n,3);
    for i in CartesianIndices(Phases)
        iz = Phases[i]
        I = CartesianIndex(i,iz)
        PhaseRatio[I] = 1.0
    end

    Hr   =  zeros(size(Phases));
    z    =  ones(size(Phases))*10e3;
    args= (z=z,)

    compute_radioactive_heat!(Hr, Mat_tup, Phases, args)
    @test  minimum(Hr)  ≈ 0.0
    @test  maximum(Hr)  ≈ 1e-6
    @test  Hr[50,50,50] ≈ 3.678794411714423e-7

    args1 = (;)
    compute_radioactive_heat!(Hr, Mat_tup1, Phases, args1)
    @test  Hr[50,50,50] ≈ 1e-6

    num_alloc = @allocated compute_radioactive_heat!(Hr, Mat_tup, Phases, args)
    @test num_alloc <= 32   # in the commandline this gives 0; while running the script not always

    compute_radioactive_heat!(Hr, Mat_tup, PhaseRatio, args)
    @test sum(Hr)≈ 0.33715177646857664

    # -----------------------

    # Shear heating -------
    Χ       = ConstantShearheating(1.0)
    @test isbits(Χ)

    # Define parameters as vectors
    τ       = [1 2 3 4]*1e6
    ε       = [1 0.1 0.1 1]
    ε_el    = [0.01 0.01 0.01 0.01]

    τ_2D       = [1 2; 3 4]*1e6
    ε_2D       = [1 0.1; 0.1 1]
    ε_el_2D    = [0.01 0.01; 0.01 0.01]

    # With elasticity
    H_s1 = compute_shearheating(Χ, τ,   ε,    ε_el   )
    H_s2 = compute_shearheating(Χ, τ_2D,ε_2D, ε_el_2D)
    @test H_s1 ≈ 5.4e6
    @test H_s2 ≈ 5.4e6

    # No elasticity
    H_s3 = compute_shearheating(Χ, τ,   ε   )
    H_s4 = compute_shearheating(Χ, τ_2D,ε_2D)
    @test H_s3 ≈ 5.5e6
    @test H_s4 ≈ 5.5e6

    # Now in non-dimensional units
    τ       = [1 2 3 4]
    ε       = [1 0.1 0.1 1]
    ε_el    = [0.01 0.01 0.01 0.01]

    τ_2D       = [1 2; 3 4]
    ε_2D       = [1 0.1; 0.1 1]
    ε_el_2D    = [0.01 0.01; 0.01 0.01]
    Χ       = nondimensionalize(Χ,CharUnits_GEO)
    info    = param_info(Χ)

    H_s1 = compute_shearheating(Χ, τ,   ε,    ε_el,  )
    H_s2 = compute_shearheating(Χ, τ_2D,ε_2D, ε_el_2D)
    H_s3 = compute_shearheating(Χ, τ,   ε)
    H_s4 = compute_shearheating(Χ, τ_2D,ε_2D)
    @test H_s1 ≈ 5.4
    @test H_s2 ≈ 5.4
    @test H_s3 ≈ 5.5
    @test H_s4 ≈ 5.5
    # -----------------------

end


text:
###
# Solve an inverse traffic problem over polynomials of degree at most d
# Optionally use a regularizer from the poly kernel
###
using JuMP
using Gurobi
using Graphs
using Roots

polyEval(coeffs, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])
polyEval(coeffs::Array{Float64, 1}, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)]) #separate for consts

bpacost(flow::Float64, capacity::Float64, freeflowtime::Float64) = freeflowtime*(1 + .15 * (flow/capacity)^4)
bpacost(flow::Float64, arc) = bpacost(flow, arc.capacity, arc.freeflowtime)
bpacost(arc::Arc) = bpacost(arc.obsflow, arc)

derivcost(flow::Float64, capacity::Float64, freeflowtime::Float64) = freeflowtime / capacity * .15 * 4 * (flow/capacity)^3
derivcost(flow::Float64, arc) = derivcost(flow, arc.capacity, arc.freeflowtime)
derivcost(arc::Arc) = derivcost(arc.obsflow, arc)

function setUpFitting(deg::Int, c, odpairs, nodes)
	m = Model(solver=GurobiSolver(OutputFlag=false))
	@defVar(m, coeffs[1:deg + 1])
	@defVar(m, Calphas[1:deg + 1])

	#build the graham matrix
	samples = linspace(0, 1, deg + 1)
	k(x,y) = (c + x*y)^deg
	K = [k(x,y) for x = samples, y=samples]
	assert(rank(K) == deg + 1)

    #Perturb a little for numerical stability
    C = chol(K + 1e-6 * eye(deg + 1))

    for i = 1:deg + 1
		@addConstraint(m, polyEval(coeffs, samples[i]) ==
						sum{C[j, i] * Calphas[j], j=1:deg+1})
	end
	@defVar(m, reg_term >= 0)
	reg_term_ = QuadExpr(Calphas[:], Calphas[:], ones(deg + 1), AffExpr() )
	addConstraint(m, reg_term >= reg_term_)

	return m, coeffs, reg_term
end

function fixCoeffs(m, fcoeffs, coeffs)
	for (fc, c) in zip(fcoeffs, coeffs[:])
		@addConstraint(m, fc == c)
	end
end

# function calcTrueCost(coeffs::Array{Float64, 1}, scaled_flows, free_flow_time, capacities)
# 	return sum([free_flow_time[a]*capacities[a]*scaled_flows[a]*polyEval(coeffs, scaled_flows[a])
# 				for a=1:length(scaled_flows)])
# end

function addResid(m, coeffs, ys, demands, arcs, scaling)
	@defVar(m, resid)
	@defVar(m, dual_cost)
	@defVar(m, primal_cost)

	@addConstraint(m, dual_cost == sum{demands[(s,t)] * (ys[(s,t), t] - ys[(s,t), s]), (s,t)=keys(demands)})
	@addConstraint(m, primal_cost == sum{a.obsflow * a.freeflowtime * polyEval(coeffs, a.obsflow/a.capacity), a=values(arcs)})

	@addConstraint(m, resid >= (dual_cost - primal_cost) / scaling )
	@addConstraint(m, resid >= (primal_cost - dual_cost) / scaling )
	return resid
end

function addIncreasingCnsts(m, coeffs, arcs; TOL=0.)
	sorted_flows = sort([a.obsflow / a.capacity for a in values(arcs)])
	@addConstraint(m, polyEval(coeffs, 0) <= polyEval(coeffs, sorted_flows[1]))
	for i = 2:length(sorted_flows)
		@addConstraint(m, polyEval(coeffs, sorted_flows[i-1]) <= polyEval(coeffs, sorted_flows[i]) + TOL)
	end
end

#equates the total cost of the network to the true total cost
function normalize(m, coeffs, tot_true_cost::Float64, arcs)
	@addConstraint(m,
		sum{a.freeflowtime * a.obsflow * polyEval(coeffs, a.obsflow / a.capacity),
			a=values(arcs)} == tot_true_cost)
end

function normalize(m, coeffs, scaled_flow::Float64, cost::Float64)
	@addConstraint(m, polyEval(coeffs, scaled_flow) == cost)
end

function normalize(m, coeffs, scaled_flows::Array{Float64, 1}, avgCost::Float64)
    @addConstraint(m, sum{polyEval(coeffs, f), f=scaled_flows} == avgCost * length(scaled_flows))
end

function addNetworkCnsts(m, coeffs, demands, arcs, numNodes)
	@defVar(m, ys[keys(demands), 1:numNodes])
	for k = keys(arcs)
		a = arcs[k]
		rhs = a.freeflowtime * polyEval(coeffs, a.obsflow/a.capacity)
		for od in keys(demands)
			@addConstraint(m, ys[od, k[2]] - ys[od, k[1]] <= rhs)
		end
	end
	return ys
end

#Uses a Frank-Wolfe algorithm to solve bpa cost for the given network.
#Takes as input the underlying graph object
#The ordering of arcs must correspond to order that the graph was constructed
function frank_wolfe(g, vArcs, demand_data, idx; TOL=1e-4, MAX_ITERS=100)
    #use the observed flows as the starting point
    flows =[a.obsflow::Float64 for a in vArcs]
    costs = [bpacost(a) for a in vArcs]
    trace = Float64[]
    for iter = 1:MAX_ITERS
        flow_dict = Dict{(Int, Int), Float64}()
        for odpair = keys(demand_data)
            #solve the shortest path problems, and update the total flow
            r = dijkstra_shortest_paths(g, costs, odpair[1] )
            currNode = odpair[2];
            while currNode != odpair[1]
                parent = r.parents[currNode]
                if ! haskey(flow_dict, (parent, currNode) )
                    flow_dict[(parent, currNode)] = demand_data[odpair][idx]
                else
                    flow_dict[(parent, currNode)] += demand_data[odpair][idx]
                end
                currNode = parent
            end
        end

        d = [get(flow_dict, (a.initNode, a.termNode), 0.)::Float64 for a in vArcs]

        #In the first iteration, just pull out the flows
        if iter == 1
        	flows = d
        	costs = [bpacost(flows[ix], a) for (ix, a) in enumerate(vArcs)]
        	continue
        end
        assert( dot(costs, d) <= dot(flows, costs) )
        d -= flows
        derivFun(alpha) = sum([bpacost(flows[ix] + alpha*d[ix], a)*d[ix] for (ix, a) in enumerate(vArcs)])
        if derivFun(0) >=0
            alpha = 0
        elseif derivFun(1) <= 0
            alpha = 1
        else
            alpha = fzero(derivFun, 0, 1)
        end
        converge_dist = alpha * norm(d) / norm(flows)
        flows += alpha * d
        push!(trace, converge_dist)
        if (iter > 1) & (converge_dist <= TOL)
            break
        else
            costs = [bpacost(flows[ix], a) for (ix, a) in enumerate(vArcs)]
        end
    end

    return trace[length(trace)], flows
end

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs


text:
module ZhuyinPinyin

export to_pinyin, to_zhuyin, Zhuyin, Pinyin, is_possible

using DataStructures

include("tries.jl")
include("pronunciations.jl")

function _replace_by_trie(source, trie)
    dest = String[]

    since_last_start = Char[]
    current_trie = trie
    for char in source
        if char in keys(current_trie.children)
            push!(since_last_start, char)
            current_trie = current_trie.children[char]
        else
            push!(dest, current_trie.is_key ? current_trie.value : join(since_last_start))
            empty!(since_last_start)

            current_trie = trie

            if char in keys(current_trie.children)
                push!(since_last_start, char)
                current_trie = current_trie.children[char]
            else
                push!(dest, string(char))
            end
        end
    end

    push!(dest, current_trie.is_key ? current_trie.value : join(since_last_start))

    return join(dest)
end

"""
    to_zhuyin(pinyin)

Convert any pinyin sequences in the source string to its zhuyin equivalent.

## Examples
```julia-repl
julia> to_zhuyin("lu4 shi1")
"ㄌㄨ4 ㄕ1"

julia> to_zhuyin("shou4 dao4 wei3 qu5")
"ㄕㄡ4 ㄉㄠ4 ㄨㄟ3 ㄑㄩ5"
```
"""
to_zhuyin(pinyin) = _replace_by_trie(pinyin, zhuyin_trie)
to_zhuyin(pinyin::Pinyin) = _replace_by_trie(pinyin.contents, zhuyin_trie)
to_zhuyin(pinyin::Zhuyin) = pinyin.contents

"""
    to_pinyin(zhuyin)

Convert any zhuyin sequences in the source string to its pinyin equivalent.

## Examples
```julia-repl
julia> to_pinyin("ㄨ ㄑㄧ ㄇㄚ ㄏㄟ")
"wu qi ma hei"

julia> to_pinyin("ㄏㄨ2 ㄕㄨㄛ1 ㄅㄚ1 ㄉㄠ4")
"hu2 shuo1 ba1 dao4"
```
"""
to_pinyin(zhuyin) = _replace_by_trie(zhuyin, pinyin_trie)
to_pinyin(zhuyin::Zhuyin) = _replace_by_trie(zhuyin.contents, pinyin_trie)
to_pinyin(zhuyin::Pinyin) = zhuyin.contents

"""
    is_possible(pronunciation)
"""
is_possible(pinyin::Pinyin) = haskey(zhuyin_trie, pinyin.contents)
is_possible(zhuyin::Zhuyin) = haskey(pinyin_trie, zhuyin.contents)

end


text:
# %j is number of days in year
time_step = Dict('s' => (Dates.Millisecond, 3), 'S' => (Dates.Second, 2),
                 'M' => (Dates.Minute, 2), 'H' => (Dates.Hour, 2),
                 'd' => (Dates.Day, 2), 'j' => (Dates.Day, 3),
                 'U' => (Dates.Month, 1), 'u' => (Dates.Month, 1), 'm' => (Dates.Month, 2),
                 'Y' => (Dates.Year, 1), 'y' => (Dates.Year, 2))
date_steps = ['s', 'S', 'M', 'H',
              'd', 'j', 'U', 'u', 'm', 'Y', 'y']
# ^^ Not using keys(time_step) because dict is not ordered. #FIXME change to ordict

function smaller_step{T <: String}(pattern::T)
    # list of date steps in increasing order
    # Loop over all the possibles step and
    # return the smallest one.
    for step in date_steps
        if step in pattern
            return step
        end
    end
    error("Date pattern $(pattern) not encountered")
end

function extract_date{T <: String}(text::T)
    # Look up for %X patterns to convert to date
    dates = collect(matchall(r"%.", text)) |> join
    replace(dates, "%", "")
end

function check_steps{T <: String}(text::T)
    # Check whether all %pattern are dates
    all(x->x in date_steps, text)
end

"""
    dates_directories(urlpattern)

Extracts the dates patterns from an url pattern
# Example
```jldoctest
julia> dates_directories("http://this.is/my/%Y/%m/%d/ofBirth.txt")
"Ymd"
```
"""
function dates_directories{T <: String}(urlpattern::T)
    directories = collect(matchall(r"%.*/", urlpattern)) |> join
    extract_date(directories)
end

"""
    create_urls(urlpattern, timeStart, timeEnd)

Generates urls directories where to look for a particular time
range.
# Example
```jldoctest
julia> dt0, dt1 = DateTime(2017, 2, 3), DateTime(2017, 2, 5)
julia> create_urls("http://this.is/my/%Y/%m/%d/ofBirth.txt", dt0, dt1)
3-element Array{String, 1}:
 "http://this.is/my/2017/02/03/"
 "http://this.is/my/2017/02/04/"
 "http://this.is/my/2017/02/05/"
```
"""
function create_urls{T <: String, DT <: Base.Dates.TimeType}(urlpattern::T, timeStart::DT, timeEnd::DT)
    #NOTE Should the steprange preceed over urlpattern?

    # Find directory step range
    dates = dates_directories(urlpattern)
    if check_steps(dates)
        step = smaller_step(dates)
    end

    # Only extract the directory path
    urlpattern_dir = matchall(r".*/", urlpattern)[1]
    function convert_url(tstep)
        # Multiple replace at once, when %j using the custom format.
        # Need to take the second element of each match (x[2]) because "%" is the first one
        replace(urlpattern_dir, Regex("%[$(join(keys(time_step)))]"),
                x -> x == "%j" ? dec(Dates.dayofyear(tstep), 3) : Dates.format(tstep, "$(x[2])" ^ time_step[x[2]][2]))
    end

    dt = timeStart:time_step[step][1](1):timeEnd
    map(convert_url, dt)

end


text:
# Unit  tests for bigfloats. Type stability of methods.

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

function bigfloats()
    @info "Test NEP-test with BigFloat"

    # Create a bigfloat NEP
    T = BigFloat
    A0 = ones(T,4,4) - Matrix{T}(1.0I, 4, 4)
    u = Vector{T}(1:4)
    v = u .- 2
    A1 = u*v'
    A2 = Matrix{T}(1.0I, 4, 4)
    A2[2,1] = T(pi)

    nep = PEP([A0,A1,A2])

    # Run augnewton once to get an accurate solution.
    v0 = ones(T, 4)
    λ0 = T(0.2)
    global itercount_star = 0
    global λiterates_star = zeros(T, 100)
    myerrmeasure = (λ,v) -> begin
        global itercount_star, λiterates_star
        itercount_star += 1
        λiterates_star[itercount_star] = λ
        return norm(compute_Mlincomb(nep, λ, v))
    end
    @info "Bigfloat precomputation"
    λstar,vstar = augnewton(T, nep; v=v0, λ=λ0, tol=eps(T)*100, errmeasure=myerrmeasure)

    @bench @testset "BigFloat comparison w $T" for T in (Float16, ComplexF16, ComplexF64)
        nep1 = PEP(Array{Array{T,2},1}(nep.A))
        global itercount = 0
        global λiterates = zeros(T, 100)
        myerrmeasure = (λ,v) -> begin
            global itercount, λiterates
            itercount += 1
            λiterates[itercount] = λ
            return abs(T(λstar - λ))
        end

        v0 = ones(T, 4)
        λ0 = T(0.2)
        λ,v = augnewton(T, nep1; v=v0, λ=λ0, tol=eps(real(T))*2, errmeasure=myerrmeasure)
        @test isa(λ, T) # Check output type
        @test abs(T(λ - λstar)) < eps(real(T))*8 # Check that we have a solution

        # Check that the kth iterate is the same as the bigfloat iteration
        k = 4
        @test abs(T(λiterates_star[k] - λiterates[k])) < eps(real(T))*10
    end
end

bigfloats()


text:
# ---
# title: 855. Exam Room
# id: problem855
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Ordered Map
# link: <https://leetcode.com/problems/exam-room/description/>
# hidden: true
# ---
#
# In an exam room, there are `N` seats in a single row, numbered `0, 1, 2, ...,
# N-1`.
#
# When a student enters the room, they must sit in the seat that maximizes the
# distance to the closest person.  If there are multiple such seats, they sit in
# the seat with the lowest number.  (Also, if no one is in the room, then the
# student sits at seat number 0.)
#
# Return a class `ExamRoom(int N)` that exposes two functions: `ExamRoom.seat()`
# returning an `int` representing what seat the student sat in, and
# `ExamRoom.leave(int p)` representing that the student in seat number `p` now
# leaves the room.  It is guaranteed that any calls to `ExamRoom.leave(p)` have
# a student sitting in seat `p`.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["ExamRoom","seat","seat","seat","seat","leave","seat"], [[10],[],[],[],[],[4],[]]
#     Output: [null,0,9,4,2,null,5]
#     **Explanation** :
#     ExamRoom(10) -> null
#     seat() -> 0, no one is in the room, then the student sits at seat number 0.
#     seat() -> 9, the student sits at the last seat number 9.
#     seat() -> 4, the student sits at the last seat number 4.
#     seat() -> 2, the student sits at the last seat number 2.
#     leave(4) -> null
#     seat() -> 5, the student sits at the last seat number 5.
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 10^9`
#   2. `ExamRoom.seat()` and `ExamRoom.leave()` will be called at most `10^4` times across all test cases.
#   3. Calls to `ExamRoom.leave(p)` are guaranteed to have a student currently sitting in seat number `p`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
quadratic_bowl_3D_lc = environment_3D(x -> transpose(x[1:2]) * x[1:2], cone = LinearizedCone)
quadratic_bowl_3D_nc = environment_3D(x -> transpose(x[1:2]) * x[1:2], cone = NonlinearCone)

circular_bowl_3D_nc = environment_3D(x -> -sqrt(2.5^2.0 - x[1]^2.0 - x[2]^2.0) + 2.5, cone = NonlinearCone)

# r = 2.0
# t = range(-r, stop = r, length = 1000)
# z = -sqrt.(r^2.0 .- t.^2.0 .- 0^2.0) .+ r
#
# plot(t, z, aspect_ratio = :equal)
# env = deepcopy(circular_bowl_3D_nc)
# env = deepcopy(flat_3D_nc)
# q = [-0.5; -0.5; env.surf([-0.5; -0.5])]
# n = [env.surf_grad(q[1:2]); -1.0]
# ns = n ./ sqrt(transpose(n) * n) #* model.r


text:
const USE_GPU = true
using ImplicitGlobalGrid
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(package=CUDA, ndims=3);
else
    @init_parallel_stencil(package=Threads, ndims=3);
end

@parallel function diffusion3D_step!(T2::Data.Array{_T}, T::Data.Array{_T}, Ci::Data.Array{_T}, lam::_T, dt::_T, _dx::_T, _dy::_T, _dz::_T) where _T <: PSNumber
    @inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
    return
end

function diffusion3D(lam::_T, c0::_T, lx::_T, ly::_T, lz::_T) where _T <: PSNumber
    # Numerics
    nx, ny, nz = 512, 512, 512;                              # Number of gridpoints in dimensions x, y and z
    nt         = 100;                                        # Number of time steps
    me, dims   = init_global_grid(nx, ny, nz);
    dx         = lx/(nx_g()-1);                              # Space step in x-dimension
    dy         = ly/(ny_g()-1);                              # Space step in y-dimension
    dz         = lz/(nz_g()-1);                              # Space step in z-dimension
    _dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;

    # Array initializations
    T   = @zeros(_T, nx, ny, nz);
    T2  = @zeros(_T, nx, ny, nz);
    Ci  = @zeros(_T, nx, ny, nz);

    # Initial conditions
    Ci .= 1/c0;                                              # 1/Heat capacity
    T  .= 1.7;
    T2 .= T;                                                 # Assign also T2 to get correct boundary conditions.

    # Time loop
    dt   = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1;          # Time step for 3D Heat diffusion
    lam, dt, _dx, _dy, _dz = _T.((lam, dt, _dx, _dy, _dz))   # Convert all scalars to the required numbertype.
    for it = 1:nt
        if (it == 11) tic(); end                             # Start measuring time.
        @parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
        update_halo!(T2);
        T, T2 = T2, T;
    end
    time_s = toc()

    # Performance
    A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(_T);               # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    t_it  = time_s/(nt-10);                                  # Execution time per iteration [s]
    T_eff = A_eff/t_it;                                      # Effective memory throughput [GB/s]
    if (me==0) println("time_s=$time_s T_eff=$T_eff"); end

    finalize_global_grid();
end

################################################################################
# CALL OF THE FUNCTION

# Physics
lam        = 1.0;                                        # Thermal conductivity
c0         = 2.0;                                        # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0;                              # Length of computational domain in dimension x, y and z

lam, c0, lx, ly, lz = Float64.((lam, c0, lx, ly, lz))    # Convert all scalars to the desired numbertype.

diffusion3D(lam, c0, lx, ly, lz)


text:
# ------------------------
# BidAsk
# ------------------------

@inline midprice(ba::BidAsk) = (ba.bid + ba.ask) / 2.0
@inline midprice(bid::Price, ask::Price) = (bid + ask) / 2.0

@inline spread(ba::BidAsk) = (ba.ask - ba.bid)
@inline spread(bid::Price, ask::Price) = (ask - bid)

@inline open_price(dir::TradeDir, ba::BidAsk) = is_long(dir) ? ba.ask : (is_short(dir) ? ba.bid : NaN)
@inline close_price(dir::TradeDir, ba::BidAsk) = is_long(dir) ? ba.bid : (is_short(dir) ? ba.ask : NaN)

# ------------------------
# TradeDir
# ------------------------

@inline is_long(dir::TradeDir) = dir === Long
@inline is_short(dir::TradeDir) = dir === Short

# ------------------------
# Position
# ------------------------

@inline is_long(pos::Position) = pos.dir === Long
@inline is_short(pos::Position) = pos.dir === Short

# size negative for shorts, thus works for both long and short
@inline pnl_net(pos::Position) = pos.size * (pos.last_price - pos.open_price)
@inline pnl_gross(pos::Position) = pos.size * (midprice(pos.last_quote) - midprice(pos.open_quote))

# size negative for shorts, thus works for both long and short
@inline return_net(pos::Position) = sign(pos.size) * (pos.last_price - pos.open_price) / pos.open_price
@inline return_gross(pos::Position) = sign(pos.size) * ((midprice(pos.last_quote) - midprice(pos.open_quote)) / midprice(pos.open_quote))

# ------------------------
# Account
# ------------------------

@inline balance_ret(acc::Account) = acc.balance / acc.initial_balance - 1.0
@inline equity_ret(acc::Account) = acc.equity / acc.initial_balance - 1.0

@inline has_open_positions(acc::Account) = length(acc.open_positions) > 0
@inline has_closed_positions(acc::Account) = length(acc.closed_positions) > 0

@inline total_return_net(acc::Account) = sum(map(return_net, acc.closed_positions))
@inline total_return_gross(acc::Account) = sum(map(return_gross, acc.closed_positions))

@inline total_pnl_net(acc::Account) = sum(map(pnl_net, acc.closed_positions))
@inline total_pnl_gross(acc::Account) = sum(map(pnl_gross, acc.closed_positions))

@inline count_winners_net(acc::Account) = count(map(x -> pnl_net(x) > 0.0, acc.closed_positions))
@inline count_winners_gross(acc::Account) = count(map(x -> pnl_gross(x) > 0.0, acc.closed_positions))

# # Dates.func(nbbo.dt) accessor shortcuts, e.g. year(nbbo), day(nbbo), hour(nbbo)
# for func in (:year, :month, :day, :hour, :minute, :second, :millisecond, :microsecond, :nanosecond)
#     name = string(func)
#     @eval begin
#         $func(ba::BidAsk)::Int64 = Dates.$func(ba.dt)
#     end
# end


text:
using Utils
using Underscores
using Chain

const INNER_DIM = 10
const OUTER_DIM = 12
const HEADER = "Tile "

struct Tile
    id::Int
    pix::Matrix{Bool}
end

function parse_file(f)
    tiles = split_blocks(f)

    re = Tile[]
    for tile in tiles
        id = 0
        pixels = falses(INNER_DIM, INNER_DIM)
        for (i,line) in enumerate(eachline(IOBuffer(tile)))
            if i == 1
                id = parse(Int, line[length(HEADER):end-1])
            else
                pixels[i-1,:] = @_ map(_ == '#', collect(line))
            end
        end
        push!(re, Tile(id, pixels))
    end
    re
end

function load()
    open(joinpath(@__DIR__, "20.txt"), "r") do f
        parse_file(f)
    end
end

tiles = load()

horizontal = rotr90 ∘ permutedims
vertical = rotl90 ∘ permutedims
const MIRRORING = [identity, horizontal, vertical]
const ROTATING = [identity, rotl90, rotr90, rot180]
all_transformations = @_ Iterators.product(MIRRORING, ROTATING) |> map(_[1] ∘ _[2], __)

function assemble(tiles)
    len = length(tiles)

    neighbors = falses(fill(len, 2)...)
    function neighborsofin(who, allowed::Set{UInt8})
        ok = typeof(allowed)()
        for (idx, neighbor) in enumerate(neighbors[:,who])
            neighbor && idx in allowed && push!(ok, idx)
        end
        ok
    end

    function connect(left, right)
        neighbors[left, right] = true
        neighbors[right, left] = true
    end

    # Find all connections.
    for left in 1:len
        l = tiles[left].pix
        for right in left+1:len
            for rot in all_transformations
                r = rot(tiles[right].pix)

                if l[INNER_DIM,:] == r[1,:]
                    connect(left, right)
                elseif l[1,:] == r[INNER_DIM,:]
                    connect(left, right)
                elseif l[:,INNER_DIM] == r[:,1]
                    connect(left, right)
                elseif l[:,1] == r[:,INNER_DIM]
                    connect(left, right)
                end
            end
        end
    end

    # Sort into categories.
    corners = Set{UInt8}()
    edges = Set{UInt8}()
    inner = Set{UInt8}()
    for i in 1:len
        ones = count(neighbors[:,i])
        to = if ones == 2
            corners
        elseif ones == 3
            edges
        else
            inner
        end
        push!(to, i)
    end

    @assert length(corners) == 4
    @assert (length(edges) / 4 + 2)^2 == len == OUTER_DIM^2

    function deduce(known, from)
        pop!(from, Iterators.only(neighborsofin(known, from)))
    end

    image = Matrix{Int}(undef, OUTER_DIM, OUTER_DIM)
    # Arbitrarily fix upper left OBDA.
    image[1,1] = pop!(corners)
    # Also the choice of these two.
    first2edges = neighborsofin(image[1,1], edges)
    image[1,2] = pop!(first2edges)
    image[2,1] = pop!(first2edges)
    delete!(edges, image[1,2])
    delete!(edges, image[2,1])

    # Deduce left and top line.
    for i in 3:OUTER_DIM-1
        image[1,i] = deduce(image[1,i-1], edges)
        image[i,1] = deduce(image[i-1,1], edges)
    end

    # The remaining corners now become clear.
    image[1,OUTER_DIM] = deduce(image[1,OUTER_DIM-1], corners)
    image[OUTER_DIM,1] = deduce(image[OUTER_DIM-1,1], corners)
    image[OUTER_DIM,OUTER_DIM] = Iterators.only(corners)

    image[OUTER_DIM,2] = deduce(image[OUTER_DIM,1], edges)
    image[2,OUTER_DIM] = deduce(image[1,OUTER_DIM], edges)
    for i in 3:OUTER_DIM-1
        image[OUTER_DIM,i] = deduce(image[OUTER_DIM,i-1], edges)
        image[i,OUTER_DIM] = deduce(image[i-1,OUTER_DIM], edges)
    end
    # All edges used.
    @assert isempty(edges)

    # Fill the inner tiles.
    for col in 2:OUTER_DIM-1, row in 2:OUTER_DIM-1
        ns = neighborsofin(image[row-1,col], inner)
        intersect!(ns, neighborsofin(image[row,col-1], inner))
        image[row,col] = pop!(inner, Iterators.only(ns))
    end
    # All inner tiles used.
    @assert isempty(inner)

    # And make a Tile matrix from the index matrix.
    reshape(map(idx -> tiles[idx], image), size(image))
end

assembled = assemble(tiles)

function orient(image::Matrix{Tile})
    pixels = Matrix{Matrix{Bool}}(undef, OUTER_DIM, OUTER_DIM)

    # The first three have to be transformed together.
    for c in all_transformations, a in all_transformations, b in all_transformations
        upper = c(image[1,1].pix)
        lower = a(image[2,1].pix)
        right = b(image[1,2].pix)
        if upper[INNER_DIM,:] == lower[1,:] && upper[:,INNER_DIM] == right[:,1]
            pixels[1,1] = upper
            pixels[2,1] = lower
            pixels[1,2] = right
            break
        end
    end

    # After that, it should be easy sailing for the rest of the row...
    for row in 3:OUTER_DIM
        prev_row = pixels[row-1,1][INNER_DIM,:]
        for t in all_transformations
            lower = t(image[row,1].pix)
            if lower[1,:] == prev_row
                pixels[row,1] = lower
                break
            end
        end
    end

    # ... and column.
    for col in 3:OUTER_DIM
        prev_col = pixels[1,col-1][:,INNER_DIM]
        for t in all_transformations
            right = t(image[1,col].pix)
            if right[:,1] == prev_col
                pixels[1,col] = right
                break
            end
        end
    end

    # Now we can always count on an upper and left neighbor.
    for col in 2:OUTER_DIM, row in 2:OUTER_DIM
        prev_row = pixels[row-1,col][INNER_DIM,:]
        prev_col = pixels[row,col-1][:,INNER_DIM]
        for t in all_transformations
            lower_right = t(image[row,col].pix)
            if lower_right[1,:] == prev_row && lower_right[:,1] == prev_col
                pixels[row,col] = lower_right
                break
            end
        end
    end

    # Lastly, pack it into a smaller image without borders.
    side = OUTER_DIM * (INNER_DIM - 2)
    re = falses(fill(side, 2)...)
    bcol = 1
    for col in 1:OUTER_DIM
        brow = 1
        for row in 1:OUTER_DIM
            re[brow:brow+INNER_DIM-3,bcol:bcol+INNER_DIM-3] = pixels[row,col][2:end-1,2:end-1]
            brow += INNER_DIM - 2
        end
        bcol += INNER_DIM - 2
    end
    re
end

oriented = orient(assembled)

const MONSTER = [
    "                  # ",
    "#    ##    ##    ###",
    " #  #  #  #  #  #   ",
]

monster = @chain MONSTER begin
    map(line -> [c == '#' for c in line], _)
    reduce(hcat, _)
    permutedims
    BitArray
end

function search(image::BitArray{2})
    rm, cm = size(monster)
    ri, ci = size(image)
    local stanced::BitArray{2}

    for t in all_transformations
        transd = t(image)
        stanced = copy(transd)

        for col in 1:ci-cm+1, row in 2:ri-1
            patch = transd[row-1:row+1,col:col+cm-1]
            if patch .& monster == monster
                # Found one, so stance it out.
                stanced[row-1:row+1,col:col+cm-1] = patch .!= monster
            end
        end

        stanced != transd && break
    end

    stanced
end

function ex1(tiles::Matrix{Tile})
    s = axes(tiles, 1)
    @chain [firstindex(s), lastindex(s)] begin
        # Twice that combination...
        fill(2)
        # ... then the combinatoric product of them...
        Iterators.product(_...)
        # ... and the multiplication of the IDs.
        prod(((r,c),) -> tiles[r,c].id, _)
    end
end

function ex2(oriented)
    oriented |> search |> count
end

println("Product of corner tiles: ", ex1(assembled))
println("Roughness of the sea: ", ex2(oriented))

using Test

@testset "Jurassic Jigsaw" begin
    @testset "results" begin
        @test ex1(assembled) == 15003787688423
        @test ex2(oriented) == 1705
    end
end

max_stars_repo_name: CmdQ/AoC


text:
"""
Return the dimension of an `Interpolation`
"""
@inline getdim(ip::Interpolation{dim}) where {dim} = dim

"""
Return the reference shape of an `Interpolation`
"""
@inline getrefshape(ip::Interpolation{dim,shape}) where {dim,shape} = shape

"""
Return the polynomial order of the `Interpolation`
"""
@inline getorder(ip::Interpolation{dim,shape,order}) where {dim,shape,order} = order

@inline gettopology(ip::Interpolation) = Dict{Int,Int}()

"""
Compute the value of the shape functions at a point ξ for a given interpolation
"""
function value(ip::Interpolation{dim}, ξ::Vec{dim,T}) where {dim,T}
    [value(ip, i, ξ) for i in 1:getnbasefunctions(ip)]
end

"""
Compute the gradients of the shape functions at a point ξ for a given interpolation
"""
function derivative(ip::Interpolation{dim}, ξ::Vec{dim,T}) where {dim,T}
    [gradient_value(ip, i, ξ) for i in 1:getnbasefunctions(ip)]
end

function value(ip::Interpolation{dim}, j::Int, ξ::Tuple) where {dim}
    value(ip,j,Vec{dim}(ξ))
end

function gradient_value(ip::Interpolation{dim}, j::Int, ξ::Tuple) where {dim}
    gradient_value(ip,j,Vec{dim}(ξ))
end

# Default to nodal interpolator for geom
function get_default_geom_interpolator(shape, ::Type{Val{1}})
    return Lagrange{1,shape,1}()
end
function get_default_geom_interpolator(shape, ::Type{Val{2}})
    return Lagrange{2,shape,1}()
end
############
# Dubiner #
###########
struct Dubiner{dim,shape,order} <: Interpolation{dim,shape,order} end

isnodal(ip::Dubiner) = false

getnbasefunctions(::Dubiner{2,RefTetrahedron,order}) where {order} = Int((order+1)*(order+2)/2)
#nvertexdofs(::Dubiner{2,RefTetrahedron,order}) = 1

#vertices(::Dubiner{2,RefTetrahedron,order}) where {order} = (1,2,3)
#faces(::Dubiner{2,RefTetrahedron,order}) where {order} = ((1,2), (2,3), (3,1))

"""
value(ip::Dubiner{2,RefTetrahedron,order}, j::Int, ξ::AbstactVector) where {order}
Compute value of dubiner basis `j` at point ξ
on the reference triangle ((0,0),(1,0),(0,1))
"""
function value(ip::Dubiner{2,RefTetrahedron,order}, j::Int, ξ::Vec{2,T}) where {order, T}
    r = ξ[1]
    s = ξ[2]
    if j == 0; return zero(T)
    elseif j == 1; return sqrt(2)*one(T)
    elseif j == 2; return 2*sqrt(3)*(2*r + s - 1)
    elseif j == 3; return 2*(3*s - 1)
    elseif j == 4; return sqrt(30)*(6*r^2 + 6*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 5; return 3*sqrt(2)*(5*s - 1)*(2*r + s - 1)
    elseif j == 6; return sqrt(6)*(10*s^2 - 8*s + 1)
    elseif j == 7; return 2*sqrt(14)*(2*r + s - 1)*(10*r^2 + 10*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 8; return 2*sqrt(10)*(7*s - 1)*(6*r^2 + 6*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 9; return 2*sqrt(6)*(21*s^2 - 12*s + 1)*(2*r + s - 1)
    elseif j == 10; return 2*sqrt(2)*(35*s^3 - 45*s^2 + 15*s - 1)
    elseif j == 11; return 3*sqrt(10)*(70*r^4 + 140*r^3*(s - 1) + 90*r^2*(s^2 - 2*s + 1) + 20*r*(s^3 - 3*s^2 + 3*s - 1) + s^4 - 4*s^3 + 6*s^2 - 4*s + 1)
    elseif j == 12; return sqrt(70)*(9*s - 1)*(2*r + s - 1)*(10*r^2 + 10*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 13; return 5*sqrt(2)*(36*s^2 - 16*s + 1)*(6*r^2 + 6*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 14; return sqrt(30)*(84*s^3 - 84*s^2 + 21*s - 1)*(2*r + s - 1)
    elseif j == 15; return sqrt(10)*(126*s^4 - 224*s^3 + 126*s^2 - 24*s + 1)
    else; return dubiner_basis(r,s,j)
    end
end

"""
gradient_value(ip::Dubiner{2,RefTetrahedron,order}, j::Int, ξ::AbstactVector) where {order}
Compute value of dubiner basis `j` derivative at point ξ
on the reference triangle ((0,0),(1,0),(0,1))
"""
function gradient_value(ip::Dubiner{2,RefTetrahedron,order}, j::Int, ξ::Vec{2,T}) where {order,T}
    r = ξ[1]
    s = ξ[2]
    if j == 0; return Vec{2}([zero(T),zero(T)])
    elseif j == 1; return Vec{2}([zero(T),zero(T)])
    elseif j == 2; return Vec{2}([4*sqrt(3), 2*sqrt(3)]*one(T))
    elseif j == 3; return Vec{2}([zero(T), one(T)*6])
    elseif j == 4; return Vec{2}([12*sqrt(30)*r + 6*sqrt(30)*(s - 1),
                            6*sqrt(30)*r + 2*sqrt(30)*s - 2*sqrt(30)])
    elseif j == 5; return Vec{2}([6*sqrt(2)*(5*s - 1),
                            30*sqrt(2)*r + 30*sqrt(2)*s - 18*sqrt(2)])
    elseif j == 6; return Vec{2}([zero(T),20*sqrt(6)*s - 8*sqrt(6)])
    elseif j == 7; return Vec{2}([24*sqrt(14)*(5*r^2 + 5*r*(s - 1) + s^2 - 2*s + 1),
                            60*sqrt(14)*r^2 + r*(48*sqrt(14)*s - 48*sqrt(14)) + 6*sqrt(14)*s^2 - 12*sqrt(14)*s + 6*sqrt(14)])
    elseif j == 8; return Vec{2}([12*sqrt(10)*(7*s - 1)*(2*r + s - 1),
                            84*sqrt(10)*r^2 + r*(168*sqrt(10)*s - 96*sqrt(10)) + 42*sqrt(10)*s^2 - 60*sqrt(10)*s + 18*sqrt(10)])
    elseif j == 9; return Vec{2}([4*sqrt(6)*(21*s^2 - 12*s + 1),
                            r*(168*sqrt(6)*s - 48*sqrt(6)) + 126*sqrt(6)*s^2 - 132*sqrt(6)*s + 26*sqrt(6)])
    elseif j == 10; return Vec{2}([zero(T),210*sqrt(2)*s^2 - 180*sqrt(2)*s + 30*sqrt(2)])
    elseif j == 11; return Vec{2}([840*sqrt(10)*r^3 + 1260*sqrt(10)*r^2*(s - 1) + 540*sqrt(10)*r*(s^2 - 2*s + 1) + 60*sqrt(10)*(s^3 - 3*s^2 + 3*s - 1),
                            420*sqrt(10)*r^3 + r^2*(540*sqrt(10)*s - 540*sqrt(10)) + r*(180*sqrt(10)*s^2 - 360*sqrt(10)*s + 180*sqrt(10)) + 12*sqrt(10)*s^3 - 36*sqrt(10)*s^2 + 36*sqrt(10)*s - 12*sqrt(10)])
    elseif j == 12; return Vec{2}([12*sqrt(70)*(9*s - 1)*(5*r^2 + 5*r*(s - 1) + s^2 - 2*s + 1),
                            180*sqrt(70)*r^3 + r^2*(540*sqrt(70)*s - 300*sqrt(70)) + r*(324*sqrt(70)*s^2 - 456*sqrt(70)*s + 132*sqrt(70)) + 36*sqrt(70)*s^3 - 84*sqrt(70)*s^2 + 60*sqrt(70)*s - 12*sqrt(70)])
    elseif j == 13; return Vec{2}([30*sqrt(2)*(36*s^2 - 16*s + 1)*(2*r + s - 1),
                            r^2*(2160*sqrt(2)*s - 480*sqrt(2)) + r*(3240*sqrt(2)*s^2 - 3120*sqrt(2)*s + 510*sqrt(2)) + 720*sqrt(2)*s^3 - 1320*sqrt(2)*s^2 + 690*sqrt(2)*s - 90*sqrt(2)])
    elseif j == 14; return Vec{2}([2*sqrt(30)*(84*s^3 - 84*s^2 + 21*s - 1),
                            r*(504*sqrt(30)*s^2 - 336*sqrt(30)*s + 42*sqrt(30)) + 336*sqrt(30)*s^3 - 504*sqrt(30)*s^2 + 210*sqrt(30)*s - 22*sqrt(30)])
    elseif j == 15; return Vec{2}([zero(T),504*sqrt(10)*s^3 - 672*sqrt(10)*s^2 + 252*sqrt(10)*s - 24*sqrt(10)])
    else; return ∇dubiner_basis(r,s,j)
    end
end

"""
    jacobi(x, p::Integer, α, β)
Evaluate the Legendre polynomial with parameters `α`, `β` of degree `p` at `x`
using the three term recursion [Karniadakis and Sherwin, Spectral/hp Element
Methods for CFD, Appendix A].
Author: H. Ranocha (see PolynomialBases.jl)
"""
function jacobi(x, p::Integer, α, β)
    T = typeof( (2+α+β)*x / 2 )
    a = one(T)
    b = ((2+α+β)*x + α - β) / 2
    if p <= 0
        return a
    elseif p == 1
        return b
    end

    for n in 2:p
        a1 = 2n*(n+α+β)*(2n-2+α+β)
        a2 = (2n-1+α+β)*(α+β)*(α-β)
        a3 = (2n-2+α+β)*(2n-1+α+β)*(2n+α+β)
        a4 = 2*(n-1+α)*(n-1+β)*(2n+α+β)
        a, b = b, ( (a2+a3*x)*b - a4*a ) / a1
    end
    b
end

"""
djacobi(x,n::Integer,α, β)
Evaluate the `n` order jacobi polynomial derivative
at point x∈[-1,1]
"""
function djacobi(x,n::Integer,α, β)
    T = typeof((2+α+β)*x/2)
    a = zero(T)
    if n <= 0;return a;end
    return (n+α+β+1)/2*jacobi(x,n-1,α+1,β+1)
end

"""
dubiner(x,y,n::Integer,m::Integer)
Compute the dubiner polynomial of degree `n`,`m` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function dubiner(x,y,n::Integer,m::Integer)
    #check domain
    @assert ((y>=0)&&(x>=0))&&(1>=x+y) "point not in domain"
    # Map to reference square
    ξ=abs(x) < eps() ? -one(typeof(x)) : 2*x/(1-y)-1
    η=2*y-1
    k=2*n+1
    #Compute Dubiner_nm(ξ, η)
    P=jacobi(ξ,n,0,0)*jacobi(η, m,k,0)*((1-η)/2)^n
    #normalize
    N=sqrt(2/((2*n+1)*(m+n+1)))
    return (2*P)/N
end

"""
∇dubiner(x,y,n::Integer,m::Integer)
Compute the gradient of dubiner polynomial of degree `(n,m)` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function ∇dubiner(x,y,n::Integer,m::Integer)
    #check domain
    @assert ((y>=0)&&(x>=0))&&(1>=x+y) "point not in domain"
    # Map to reference square
    ξ=abs(x) < eps() ? -one(typeof(x)) : 2*x/(1-y)-1
    k=2*n+1
    η=2*y-1
    #Compute ∇Dubiner_nm(ξ, η)
    Dφn = djacobi(ξ,n,0,0)
    φn = jacobi(ξ,n,0,0)
    φm = jacobi(η,m,k,0)
    Dφm = djacobi(η,m,k,0)
    Px=2/(1-η)*Dφn*φm*((1-η)/2)^n
    N=sqrt(2/((2*n+1)*(m+n+1)))
    Py=(2*x/(1-y)^2*Dφn*((1-η)/2)^n-n*((1-η)/2)^(n-1)*φn)*φm + 2*φn*((1-η)/2)^n*Dφm
    return [(4*Px)/N,(2*Py)/N]
end

"""
dubiner_basis(x,y,j::Integer)
Evaluate the dubiner basis `j` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function dubiner_basis(x,y,j::Integer)
    #Compute degrees
    t=-3/2+(1/2)*sqrt(1+8*j)
    n=((ceil(t)+1)*(ceil(t)+2))/2-j
    m=ceil(t)-n
    #Compute Dubiner_nm(ξ, η)
    return dubiner(x,y,Int(n),Int(m))
end

"""
∇dubiner_basis(x,y,j::Integer)
Compute the gradient of dubiner basis `j` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function ∇dubiner_basis(x,y,j::Integer)
    #Compute degrees
    t=-3/2+(1/2)*sqrt(1+8*j)
    n=Int(((ceil(t)+1)*(ceil(t)+2))/2-j)
    m=Int(ceil(t)-n)
    #Compute ∇Dubiner_nm(ξ, η)
    return Vec{2}(∇dubiner(x,y,n,m))
end

####################
# Lagrange
####################
struct Lagrange{dim,shape,order} <: Interpolation{dim,shape,order}
    nodal_base_coefs::Matrix{Float64}
    topology::Dict{Int,Int}
end

isnodal(ip::Lagrange) = true

function getdefaultdualbasis(dim::Int,shape,order::Int)
    if dim == 1
        return Legendre{1,shape,order}()
    elseif dim == 2
        return Dubiner{2,shape,order}()
    else
        throw("Not dual basis available for dimension $dim")
    end
end

function defaultDualBasis(dim::Int,shape,order::Int)
    if dim == 1
        return Legendre{1,shape,order}()
    elseif dim == 2
        return Dubiner{2,shape,order}()
    else
        throw("Not dual basis available for dimension $dim")
    end
end

function Lagrange{dim,shape,order}() where {dim,shape, order}
    nodal_points, topology = get_nodal_points(shape, Val{dim}, order)
    nodals=[x->x(point) for point in nodal_points]
    ip_prime = getdefaultdualbasis(dim,shape,order)
    nbasefuncs = getnbasefunctions(ip_prime)
    prime_base = [x->value(ip_prime, j, x) for j in 1:nbasefuncs]
    V = reshape([nodals[i](prime_base[j]) for j = 1:nbasefuncs for i=1:nbasefuncs],(nbasefuncs,nbasefuncs))
    nodal_base_coefs = inv(V)
    Lagrange{dim, shape, order}(nodal_base_coefs, topology)
end

@inline getnbasefunctions(::Lagrange{1,RefTetrahedron,order}) where {order} = order + 1
@inline getnbasefunctions(::Lagrange{2,RefTetrahedron,order}) where {order} = Int((order+1)*(order+2)/2)
@inline gettopology(ip::Lagrange{dim,RefTetrahedron,order}) where {dim,order} = ip.topology

"""
return interpolator of the same type with dim = dim -1
"""
@inline getlowerdiminterpol(ip::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim-1,shape,order}()

"""
value(ip::Lagrange{2,RefTetrahedron,order}, j::Int, ξ::AbstactVector) where {order}
Compute value of Lagrange basis `j` at point ξ
on the reference triangle ((0,0),(1,0),(0,1))
"""
function value(ip::Lagrange{2,RefTetrahedron,order}, k::Int, ξ::Vec{2,T}) where {order, T}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    n = getnbasefunctions(ip)
    dot(ip.nodal_base_coefs[:,k], value(Dubiner{2,RefTetrahedron,order}(), ξ))
end

function value(ip::Lagrange{1,RefTetrahedron,order}, k::Int, ξ::Vec{1,T}) where {order, T}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    n = getnbasefunctions(ip)
    dot(ip.nodal_base_coefs[:,k], value(Legendre{1,RefTetrahedron,order}(), ξ))
end

"""
gradient_value(ip::Lagrange{2,RefTetrahedron,order}, j::Int, ξ::AbstactVector) where {order}
Compute value of Lagrange basis `j` derivative at point ξ
on the reference triangle ((0,0),(1,0),(0,1))
"""
function gradient_value(ip::Lagrange{dim,RefTetrahedron,order}, k::Int, ξ::Vec{dim,T}) where {dim,order,T}
    if k >getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    gradient(ξ -> value(ip, k, ξ), ξ)
end

function _get_nodal_transformation_matrix(fe::Lagrange{dim,shape,order}) where {dim,shape,order}
    # Matrix to get spacial coordinates
    nodal_points, topology = get_nodal_points(shape, Val{dim}, order)
    T = eltype(nodal_points[1])
    geom_interpol = get_default_geom_interpolator(shape, Val{dim})
    qrs = QuadratureRule{dim,shape,T}(fill(T(NaN), length(nodal_points)), nodal_points) # weights will not be used
    n_qpoints = length(getweights(qrs))
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints)
    for (qp, ξ) in enumerate(qrs.points)
        for i in 1:n_geom_basefuncs
            M[i, qp] = value(geom_interpol, i, ξ)
        end
    end
    M
end

function spatial_nodal_coordinate(fe::Lagrange, M::Matrix{T}, n_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_base_funcs = size(M,1)
    @assert length(x) == n_base_funcs
    vec = zero(Vec{dim,T})
    @inbounds for i in 1:n_base_funcs
        vec += M[i, n_point] * x[i]
    end
    return vec
end

####################
# Legendre
####################
struct Legendre{dim,shape,order} <: Interpolation{dim,shape,order} end
isnodal(Legendre) = false

@inline getnbasefunctions(::Legendre{1,RefTetrahedron,order}) where {order} = order + 1

"""
value(ip::Legendre{1,RefTetrahedron,order}, j::Int, ξ::AbstactVector) where {order}
Compute value of Legendre basis `j` at point ξ
on the reference line (0,1)
"""
function value(ip::Legendre{1,RefTetrahedron,order}, k::Int, ξ::Vec{1,T}) where {order, T}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    return sqrt((2*(k-1)+1))*jacobi(2*ξ[1]-1,k-1,0.0,0.0)
end

"""
gradient_value(ip::Legendre{1,RefTetrahedron,order}, j::Int, ξ::AbstactVector) where {order}
Compute value of Legendre basis `j` derivative at point ξ
on the reference line (0,1)
"""
function gradient_value(ip::Legendre{1,RefTetrahedron,order}, k::Int, ξ::Vec{1,T}) where {order, T}
    if k >getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    gradient(ξ -> value(ip, k, ξ), ξ)
end


text:
module AtomicGraphNets

export AGNConv, AGNPool#, AGNConvDEQ
include("layers.jl")
using .Layers: AGNConv, AGNPool

include("models.jl")
export build_CGCNN, build_SGCNN

end


text:
""" AST and pretty printer for TikZ.

This module does not provide bindings to the TikZ LaTeX package. For that, see
the TikzPictures.jl package: https://github.com/sisl/TikzPictures.jl

The AST is large but still incomplete! It supports:

- Nodes (`\\node`) and edges (`\\draw`)
- Nodes along edges (`\\draw ... node ...`)
- Graphs (`\\graph`)
- Matrices (`\\matrix`)
- Scopes and nested pictures

The AST is adapted from the (also incomplete) BNF grammar for TikZ in
[TikZit](http://tikzit.sourceforge.net/manual.html).
"""
module TikZ
export Expression, PathExpression, Statement, GraphStatement, Property,
  Document, Picture, Scope, Coordinate, Node, NodeCoordinate, Edge, EdgeNode,
  PathOperation, Graph, GraphScope, GraphNode, GraphEdge, MatrixNode, pprint

using Compat
using AutoHashEquals

# AST
#####

abstract type Expression end
abstract type PathExpression <: Expression end
abstract type Statement <: Expression end
abstract type GraphStatement <: Expression end

@auto_hash_equals struct Property <: Expression
  key::String
  value::Union{String,Vector{Property},Nothing}

  Property(key::String, value=nothing) = new(key, value)
end

as_property(prop::Property) = prop
as_property(pair::Pair) = as_property(first(pair), last(pair))
as_property(key::String, value::Union{String,Nothing}=nothing) = Property(key, value)
as_property(key::String, values::AbstractVector) = Property(key, as_properties(values))

as_properties(props::Vector{Property}) = props
as_properties(props) = Property[ as_property(prop) for prop in props ]

@auto_hash_equals struct Picture <: Expression
  stmts::Vector{Statement}
  props::Vector{Property}

  Picture(stmts::Vararg{Statement}; props=Property[]) = new([stmts...], props)
end

@auto_hash_equals struct Document <: Expression
  picture::Picture
  libraries::Vector{String}
  packages::Vector{String}

  Document(picture::Picture; libraries=String[], packages=String[]) =
    new(picture, libraries, packages)
end

@auto_hash_equals struct Scope <: Statement
  stmts::Vector{Statement}
  props::Vector{Property}

  Scope(stmts::Vararg{Statement}; props=Property[]) = new([stmts...], props)
end

@auto_hash_equals struct Coordinate <: PathExpression
  x::String
  y::String

  Coordinate(x::String, y::String) = new(x, y)
  Coordinate(x::Number, y::Number) = new(string(x), string(y))
end

@auto_hash_equals struct Node <: Statement
  name::Union{String,Nothing}
  props::Vector{Property}
  coord::Union{Coordinate,Nothing}
  content::String

  Node(; name=nothing, props=Property[], coord=nothing, content="") =
    new(name, props, coord, content)
end
Node(name::Union{String,Nothing}; kw...) = Node(; name=name, kw...)

@auto_hash_equals struct NodeCoordinate <: PathExpression
  name::String
end

@auto_hash_equals struct EdgeNode <: PathExpression
  props::Vector{Property}
  content::Union{String,Nothing}

  EdgeNode(; props=Property[], content=nothing) = new(props, content)
end

@auto_hash_equals struct Edge <: Statement
  exprs::Vector{PathExpression}
  props::Vector{Property}

  Edge(args...; props=Property[]) =
    new([arg isa String ? NodeCoordinate(arg) : arg for arg in args], props)
end

@auto_hash_equals struct PathOperation <: PathExpression
  op::String
  props::Vector{Property}

  PathOperation(op::String; props=Property[]) = new(op, props)
end

@auto_hash_equals struct Graph <: Statement
  stmts::Vector{GraphStatement}
  props::Vector{Property}

  Graph(stmts::Vararg{GraphStatement}; props=Property[]) =
    new([stmts...], props)
end

@auto_hash_equals struct GraphScope <: GraphStatement
  stmts::Vector{GraphStatement}
  props::Vector{Property}

  GraphScope(stmts::Vararg{GraphStatement}; props=Property[]) =
    new([stmts...], props)
end

@auto_hash_equals struct GraphNode <: GraphStatement
  name::String
  props::Vector{Property}
  content::Union{String,Nothing}

  GraphNode(name::String; props=Property[], content=nothing) =
    new(name, props, content)
end

@auto_hash_equals struct GraphEdge <: GraphStatement
  src::String
  tgt::String
  props::Vector{Property}

  GraphEdge(src::String, tgt::String; props=Property[]) =
    new(src, tgt, props)
end

@auto_hash_equals struct MatrixNode <: Statement
  stmts::Matrix{Vector{Statement}}
  props::Vector{Property}

  MatrixNode(stmts::Matrix; props=Property[]) = new(stmts, props)
end

# Pretty-print
##############

""" Pretty-print the TikZ expression.
"""
pprint(expr::Expression; kw...) = pprint(stdout, expr; kw...)
pprint(io::IO, expr::Expression; kw...) = pprint(io, expr, 0; kw...)

function pprint(io::IO, doc::Document, n::Int)
  for package in doc.packages
    indent(io, n)
    println(io, "\\usepackage{$package}")
  end
  for library in doc.libraries
    indent(io, n)
    println(io, "\\usetikzlibrary{$library}")
  end
  pprint(io, doc.picture, n)
end

function pprint(io::IO, pic::Picture, n::Int)
  indent(io, n)
  print(io, "\\begin{tikzpicture}")
  pprint(io, pic.props)
  println(io)
  for stmt in pic.stmts
    pprint(io, stmt, n+2)
    println(io)
  end
  indent(io, n)
  print(io, "\\end{tikzpicture}")
end

function pprint(io::IO, scope::Scope, n::Int)
  indent(io, n)
  print(io, "\\begin{scope}")
  pprint(io, scope.props)
  println(io)
  for stmt in scope.stmts
    pprint(io, stmt, n+2)
    println(io)
  end
  indent(io, n)
  print(io, "\\end{scope}")
end

function pprint(io::IO, node::Node, n::Int)
  indent(io, n)
  print(io, "\\node")
  pprint(io, node.props)
  if !isnothing(node.name)
    print(io, " ($(node.name))")
  end
  if !isnothing(node.coord)
    print(io, " at ")
    pprint(io, node.coord)
  end
  print(io, " {$(node.content)};")
end

function pprint(io::IO, edge::Edge, n::Int)
  indent(io, n)
  print(io, "\\path")
  pprint(io, edge.props)
  for expr in edge.exprs
    print(io, " ")
    pprint(io, expr)
  end
  print(io, ";")
end

function pprint(io::IO, node::EdgeNode, n::Int)
  print(io, "node")
  pprint(io, node.props)
  if !isnothing(node.content)
    print(io, " {$(node.content)}")
  end
end

function pprint(io::IO, graph::Graph, n::Int)
  indent(io, n)
  print(io, "\\graph")
  pprint(io, graph.props)
  println(io, "{")
  for stmt in graph.stmts
    pprint(io, stmt, n+2)
    println(io)
  end
  indent(io, n)
  print(io, "};")
end

function pprint(io::IO, scope::GraphScope, n::Int)
  indent(io, n)
  print(io, "{")
  pprint(io, scope.props)
  println(io)
  for stmt in scope.stmts
    pprint(io, stmt, n+2)
    println(io)
  end
  indent(io, n)
  print(io, "};")
end

function pprint(io::IO, node::GraphNode, n::Int)
  indent(io, n)
  print(io, node.name)
  if !isnothing(node.content)
    print(io, "/\"$(node.content)\"")
  end
  if !isempty(node.props)
    print(io, " ")
    pprint(io, node.props)
  end
  print(io, ";")
end

function pprint(io::IO, node::GraphEdge, n::Int)
  indent(io ,n)
  print(io, "$(node.src) ->")
  pprint(io, node.props)
  print(io, " $(node.tgt);")
end

function pprint(io::IO, matrix::MatrixNode, ind::Int)
  indent(io, ind)
  print(io, "\\matrix")
  pprint(io, matrix.props)
  println(io, "{")
  m,n = size(matrix.stmts)
  for i = 1:m
    for j = 1:n
      p = length(matrix.stmts[i,j])
      for k = 1:p
        pprint(io, matrix.stmts[i,j][k], ind+2)
        if (k < p) println(io) end
      end
      if (j < n) println(io, " &") end
    end
    println(io, " \\\\")
  end
  indent(io, ind)
  print(io, "};")
end

function pprint(io::IO, coord::Coordinate, n::Int)
  print(io, "($(coord.x),$(coord.y))")
end

function pprint(io::IO, node::NodeCoordinate, n::Int)
  print(io, "($(node.name))")
end

function pprint(io::IO, op::PathOperation, n::Int)
  print(io, op.op)
  pprint(io, op.props)
end

function pprint(io::IO, prop::Property, n::Int)
  print(io, prop.key)
  if !isnothing(prop.value)
    print(io, "=")
    if prop.value isa Vector{Property}
      pprint(io, prop.value, nested=true)
    else
      print(io, prop.value)
    end
  end
end

function pprint(io::IO, props::Vector{Property}, n::Int=0; nested::Bool=false)
  if !isempty(props)
    print(io, nested ? "{" : "[")
    for (i, prop) in enumerate(props)
      pprint(io, prop)
      if (i < length(props)) print(io, ",") end
    end
    print(io, nested ? "}" : "]")
  end
end

indent(io::IO, n::Int) = print(io, " "^n)

end


text:
# This is loaded conditionally via @require if PyPlot is loaded

using .PyPlot: plt, PyCall

export pgdsgui

get_bystr(@nospecialize(by)) = by === inclusive ? "Inclusive" :
                               by === exclusive ? "Exclusive" : error("unknown ", by)

"""
    methodref, ax = pgdsgui(tinf::InferenceTimingNode; consts::Bool=true, by=inclusive)
    methodref     = pgdsgui(ax, tinf::InferenceTimingNode; kwargs...)

Create a scatter plot comparing:
    - (vertical axis) the inference time for all instances of each Method, as captured by `tinf`;
    - (horizontal axis) the run time cost, as estimated by capturing a `@profile` before calling this function.

Each dot corresponds to a single method. The face color encodes the number of times that method was inferred,
and the edge color corresponds to the fraction of the runtime spent on runtime dispatch (black is 0%, bright red is 100%).
Clicking on a dot prints the method (or location, if inlined) to the REPL, and sets `methodref[]` to
that method.

`ax` is the pyplot axis of the scatterplot.

!!! compat
    `pgdsgui` depends on PyPlot via the Requires.jl package. You must load both SnoopCompile and PyPlot for this function to be defined.
"""
function pgdsgui(ax::PyCall.PyObject, ridata; bystr, consts, markersz=25, linewidth=0.5, t0 = 0.001, interactive::Bool=true, kwargs...)
    methodref = Ref{Union{Method,MethodLoc}}()
    function onclick(event)
        xc, yc = event.xdata, event.ydata
        idx = argmin((log.(rts .+ t0) .- log(xc)).^2 + (log.(its .+ t0) .- log(yc)).^2)
        m = meths[idx]
        methodref[] = m
        println(m, " ($(nspecs[idx]) specializations)")
    end

    meths, rts, its, nspecs, ecols = Union{Method,MethodLoc}[], Float64[], Float64[], Int[], Tuple{Float64,Float64,Float64}[]
    for (m, d) in ridata  # (rt, trtd, it, nspec)
        push!(meths, m)
        push!(rts, d.trun)
        push!(its, d.tinf)
        push!(nspecs, d.nspec)
        push!(ecols, (d.trun > 0 ? d.trtd/d.trun : 0.0, 0.0, 0.0))
    end
    sp = sortperm(nspecs)
    meths, rts, its, nspecs, ecols = meths[sp], rts[sp], its[sp], nspecs[sp], ecols[sp]

    # Add t0 to each entry to handle times of zero in the log-log plot
    smap = ax.scatter(rts .+ t0, its .+ t0, markersz, nspecs; norm=plt.matplotlib.colors.LogNorm(), edgecolors=ecols, linewidths=linewidth, kwargs...)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_xlabel("Run time (self) + $t0 (s)")
    ax.set_ylabel("$bystr inference time + $t0 (s)")
    # ax.set_aspect("equal")
    ax.axis("square")
    constsmode = consts ? "incl." : "excl."
    plt.colorbar(smap, label = "# specializations ($constsmode consts)", ax=ax)
    if interactive
        ax.get_figure().canvas.mpl_connect("button_press_event", onclick)
    end
    return methodref
end

function pgdsgui(ridata; kwargs...)
    fig, ax = plt.subplots()
    return pgdsgui(ax, ridata; kwargs...), ax
end

function pgdsgui(ax::PyCall.PyObject, tinf::InferenceTimingNode; consts::Bool=true, by=inclusive, kwargs...)
    pgdsgui(ax, runtime_inferencetime(tinf; consts, by); bystr=get_bystr(by), consts, kwargs...)
end
function pgdsgui(tinf::InferenceTimingNode; consts::Bool=true, by=inclusive, kwargs...)
    pgdsgui(runtime_inferencetime(tinf; consts, by); bystr=get_bystr(by), consts, kwargs...)
end

@deprecate specialization_plot pgdsgui


text:
URLs= ["rp-f04972.local", "rp-f044d6.local"]
mkpath("images")

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
# Code for video 3_1:

fname = "iris.txt"
open(fname) do file
    for line in eachline(file)
        # process line
    end
end

type Iris
    sepal_length::Float64
    sepal_width::Float64
    petal_length::Float64
    petal_width::Float64
    species::ASCIIString
end

iris1 = Iris(5.1, 3.5, 1.4, 0.2, "I. setosa")

fname = "example.dat"
open(fname, "w") do filemode(file)
    write(file, string(iris1))
    write(file, "\n")
    # returns 33 (bytes written)
    println(file, "*** End of data ***")
end

# reading delimited files:
fname = "iris.csv"
data = readcsv(fname)
# 150x5 Array{Any,2}:
#  5.1  3.5  1.4  0.2  "I. setosa"
#  4.9  3.0  1.4  0.2  "I. setosa"
#  4.7  3.2  1.3  0.2  "I. setosa"
#  4.6  3.1  1.5  0.2  "I. setosa"
#  5.0  3.6  1.4  0.2  "I. setosa"
#  5.4  3.9  1.7  0.4  "I. setosa"
#  4.6  3.4  1.4  0.3  "I. setosa"
#  5.0  3.4  1.5  0.2  "I. setosa"
#  4.4  2.9  1.4  0.2  "I. setosa"
#  4.9  3.1  1.5  0.1  "I. setosa"
#  5.4  3.7  1.5  0.2  "I. setosa"
#  4.8  3.4  1.6  0.2  "I. setosa"
#  ⋮
#  6.0  3.0  4.8  1.8  "I. virginica"
#  6.9  3.1  5.4  2.1  "I. virginica"
#  6.7  3.1  5.6  2.4  "I. virginica"
#  6.9  3.1  5.1  2.3  "I. virginica"
#  5.8  2.7  5.1  1.9  "I. virginica"
#  6.8  3.2  5.9  2.3  "I. virginica"
#  6.7  3.3  5.7  2.5  "I. virginica"
#  6.7  3.0  5.2  2.3  "I. virginica"
#  6.3  2.5  5.0  1.9  "I. virginica"
#  6.5  3.0  5.2  2.0  "I. virginica"
#  6.2  3.4  5.4  2.3  "I. virginica"
#  5.9  3.0  5.1  1.8  "I. virginica"

data = readdlm(fname, ',')
# data = readdlm(fname, ';', Float64, '\n', header=true)

fname = "iris.dlm"
writedlm(fname, data, ';')

# HDF5 format:
# https://en.wikipedia.org/wiki/Hierarchical_Data_Format
# https://github.com/timholy/HDF5.jl
Pkg.add("HDF5")
using HDF5

h5open("iris.h5", "w") do file
    write(file, "iris1", string(iris1))
    # alternatively: @write file iris1
end

data = h5read("iris.h5", "iris1") # "Iris(5.1,3.5,1.4,0.2,\"I. setosa\")"
# alternatively:
data = h5open("iris.h5", "r") do file
    read(file, "iris1")
end

file = h5open("iris.h5", "w")
g = g_create(file, "iris") # create a group
g["dset1"] = 150           # create a scalar dataset inside the group, equivalent: write(g, "dset1", 150)
attrs(g)["Description"] = "This group contains the number of iris measurements" # an attribute
close(file)

file = h5open("iris.h5", "r") # HDF5 data file: iris.h5
g = file["iris"] # HDF5 group: /iris (file: iris.h5)
n = read(g, "dset1") # 150 - read dataset dset1 from group g

dump(file)
# HDF5.HDF5File len 1
#   iris: HDF5.HDF5Group len 1
#     dset1: HDF5Dataset () : 150

# JLD format:
# https://github.com/JuliaLang/JLD.jl
Pkg.add("JLD")
using JLD

species = "I. virginica"
save("d:\\temp\\iris.jld", "species", species) # on Windows
# save("/temp/iris.jld", "species", species) # on Linux
new = load("d:\\temp\\iris.jld")
# Dict{ByteString,Any} with 1 entry:
#   "species" => "I. virginica"
new = load("d:\\temp\\iris.jld", "species") # "I. virginica"

jldopen("d:\\temp\\iris.jld", "w") do file
    write(file, "length", 5.1)
end

data = jldopen("d:\\temp\\iris.jld", "r") do file
    read(file, "length")
end # 5.1

using HDF5
file = jldopen("d:\\temp\\iris.jld", "w")
g = g_create(file, "iris") # create a group
g["dset1"] = 150           # create a scalar dataset inside the group, equivalent: write(g, "dset1", 150)
close(file)

new = load("d:\\temp\\iris.jld")
# Dict{ByteString,Any} with 1 entry:
#   "iris" => Dict{ByteString,Any}("dset1"=>150)

# JSON
Pkg.add("JSON")
using JSON

fname = "iris.json"
data = readall(fname)
data = replace(data, "\r\n", "")
data = replace(data, " ", "")
# "{\"species\":\"I.setosa\",\"sepal_length\":5.1,\"sepal_width\":3.5,\"petal_lengt
# h\":1.4,\"petal_width\":0.2}{\"species\":\"I.setosa\",\"sepal_length\":4.9,\"sepa
# l_width\":3.0,\"petal_length\":1.4,\"petal_width\":0.2}{\"species\":\"I.virginica
# \",\"sepal_length\":5.9,\"sepal_width\":3.0,\"petal_length\":5.1,\"petal_width\":
# 1.8}"
data = replace(data, "}{", "}-{") # to split the string in an array pn '-'
data = split(data, "-")
for rec in data
    iris = JSON.parse(rec)
    @printf "Species: %s, sepal_length: %1.1f\n" iris["species"] iris["sepal_length"]
end
# Species: I.setosa, sepal_length: 5.1
# Species: I.setosa, sepal_length: 4.9
# Species: I.virginica, sepal_length: 5.9

d = Dict(:species=>"I.setosa",:sepal_length=>5.1,:sepal_width=>3.5,:petal_
length=>1.4,:petal_width=>0.2)
# Dict{Symbol,Any} with 5 entries:
#   :sepal_width  => 3.5
#   :sepal_length => 5.1
#   :petal_length => 1.4
#   :species      => "I.setosa"
#   :petal_width  => 0.2

julia> JSON.json(d)
# "{\"sepal_width\":3.5,\"sepal_length\":5.1,\"petal_length\":1.4,\"species\":\"I.s
# etosa\",\"petal_width\":0.2}"

download("https://www.googleapis.com/books/v1/volumes?q=isbn:1408855895")
# "C:\\Users\\CVO\\AppData\\Local\\Temp\\jul2BE0.tmp"
download("https://www.googleapis.com/books/v1/volumes?q=isbn:1408855895", "book.json")
#"book.json"

max_stars_repo_path: Section 3/3_1.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 3/3_1.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 3/3_1.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
@with_kw mutable struct Title
    text::String = ""
    # ---
    textstyle::TextStyle = TextStyle()
    # ---
    dist::Option{Float64} = ∅  # distance labels - title
end

@with_kw mutable struct TicksLabels
    names::Vector{String} = String[]
    # ---
    textstyle::TextStyle = TextStyle()
    # ---
    angle ::Option{Float64} = ∅ # rotation of labels
    format::Option{String}  = ∅ # format of the ticks labels
    shift ::Option{Float64} = ∅ # move labels to left/right
    dist  ::Option{Float64} = ∅ # ⟂ distance to spine
    # --- toggle-able
    off   ::Bool = false # whether to suppress the labels
end

@with_kw mutable struct Ticks
    labels   ::TicksLabels = TicksLabels() # their label
    linestyle::LineStyle   = LineStyle()   # how the ticks marks look
    # ---
    places   ::Vector{Float64} = Float64[] # where the ticks are
    length   ::Option{Float64} = ∅         # how long the ticks spine (negative for outside)
    # --- toggle-able
    symticks ::Bool = false # draws ticks on 2 sides of
    off      ::Bool = false # whether to suppress them
    grid     ::Bool = false # ticks increased to mirrorred axis
end


text:
module XAM

export
    SAM,
    BAM

"""
    flag(record::Union{SAM.Record, BAM.Record})::UInt16

Get the bitwise flags of `record`. The returned value is a `UInt16` of each flag
being OR'd together. The possible flags are:

    0x0001 template having multiple segments in sequencing
    0x0002 each segment properly aligned according to the aligner
    0x0004 segment unmapped
    0x0008 next segment in the template unmapped
    0x0010 SEQ being reverse complemented
    0x0020 SEQ of the next segment in the template being reverse complemented
    0x0040 the first segment in the template
    0x0080 the last segment in the template
    0x0100 secondary alignment
    0x0200 not passing filters, such as platform/vendor quality controls
    0x0400 PCR or optical duplicate
    0x0800 supplementary alignment
"""
function flag end

include("sam/sam.jl")
include("bam/bam.jl")

using .SAM
using .BAM

end # module


text:
"""
    DualizedNLPBlockModel{T, S} <: AbstractNLPModel{T, S}

A data type to store dualized subproblems.
"""
mutable struct DualizedNLPBlockModel{T,S} <: AbstractNLPModel{T,S}
    meta::NLPModelMeta{T,S}
    counters::Counters
    problem_block::AbstractNLPModel
    λ::AbstractVector # dual variables
    A::AbstractArray # linking matrices
end

"""
    DualizedNLPBlockModel(
        nlp::AbstractNLPModel,
        λ::AbstractVector,
        A::AbstractArray
    )

Modifies a subproblem by dualizing the linking constraints.

# Arguments
- `nlp::AbstractNLPModel`: the subproblem
- `λ::AbstractVector`: vector of dual variables
- `A::AbstractArray`: block linking matrix corresponding to the subproblem `nlp`
"""
function DualizedNLPBlockModel(nlp::AbstractNLPModel, λ::AbstractVector, A::AbstractArray)
    meta = nlp.meta
    return DualizedNLPBlockModel(meta, Counters(), nlp, λ, A)
end

"""
    update_dual!(
      nlp::DualizedNLPBlockModel,
      λ::AbstractVector,
    )
Updates the dual solution in-place for the dualized nlp block `nlp`.

# Arguments
- `nlp::DualizedNLPBlockModel`: the subproblem
- `λ::AbstractVector`: vector of dual variables
"""
function update_dual!(nlp::DualizedNLPBlockModel, λ::AbstractVector)
    nlp.λ .= λ
end

function NLPModels.obj(nlp::DualizedNLPBlockModel, x::AbstractVector)
    n = nlp.problem_block.meta.nvar
    return obj(nlp.problem_block, x) + dot(nlp.λ, nlp.A, x)
end

function NLPModels.grad!(nlp::DualizedNLPBlockModel, x::AbstractVector, g::AbstractVector)
    grad!(nlp.problem_block, x, g)
    mul!(g, nlp.A', nlp.λ, 1, 1)
    return g
end

function NLPModels.hess_structure!(
    nlp::DualizedNLPBlockModel,
    rows::AbstractVector{T},
    cols::AbstractVector{T},
) where {T}
    return hess_structure!(nlp.problem_block, rows, cols)
end

function NLPModels.hess_coord!(
    nlp::DualizedNLPBlockModel,
    x::AbstractVector{T},
    vals::AbstractVector{T};
    obj_weight = one(T),
) where {T}
    return hess_coord!(nlp.problem_block, x, vals, obj_weight = obj_weight)
end

function NLPModels.hess_coord!(
    nlp::DualizedNLPBlockModel,
    x::AbstractVector{T},
    y::AbstractVector{T},
    vals::AbstractVector{T};
    obj_weight = one(T),
) where {T}
    return hess_coord!(nlp.problem_block, x, y, vals, obj_weight = obj_weight)
end

function NLPModels.cons!(nlp::DualizedNLPBlockModel, x::AbstractVector, cx::AbstractVector)
    return cons!(nlp.problem_block, x, cx)
end

function NLPModels.jac_structure!(
    nlp::DualizedNLPBlockModel,
    rows::AbstractVector{T},
    cols::AbstractVector{T},
) where {T}
    return jac_structure!(nlp.problem_block, rows, cols)
end

function NLPModels.jac_coord!(
    nlp::DualizedNLPBlockModel,
    x::AbstractVector,
    vals::AbstractVector,
)
    return jac_coord!(nlp.problem_block, x, vals)
end


text:
### A Pluto.jl notebook ###
# v0.19.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ a68432ac-b7eb-4985-8232-67bde7bcfd0a
# ╠═╡ show_logs = false
begin
	import Pkg
    Pkg.activate(mktempdir())	#Pkg.activate("/home/jeff/Jules/PNToyBox")
	Pkg.add(path="/home/jeff/Jules/PNToyBox/src/PNML")
	Pkg.add("PlutoUI")
	Pkg.add("LabelledArrays")
	Pkg.add("OrdinaryDiffEq")
	Pkg.add("EzXML")
	Pkg.add("Plots")
	Pkg.add("Petri")

	using PlutoUI
    using PNML, LabelledArrays, OrdinaryDiffEq
	using EzXML
    using Plots
	using Petri
end

# ╔═╡ e6e05790-dd20-11ec-39b0-eb0589d76e44
md"""
# PNML Lotka-Volterra Notebook
"""

# ╔═╡ 0143155c-a570-4a03-9cab-5704b7e49968

# ╔═╡ cff98b32-c6d1-4c4d-beba-ab46f9a54172
# ╠═╡ show_logs = false
TableOfContents(title="TOC", aside=true)

# ╔═╡ 97ed66fd-e8aa-4291-81f2-c4981c521862
# ╠═╡ show_logs = false
versioninfo()

# ╔═╡ 9420e502-ed69-443b-944b-9339f2e699b1
# ╠═╡ show_logs = false
Pkg.status()

# ╔═╡ 5db3c534-07dc-4a9d-8782-41423da2ccfe
md"""
## Scratchpad
"""

# ╔═╡ a995864c-937e-4a0d-b99a-31fd134e23dc
let
	d = PNML.PnmlDict(:labels => PNML.PnmlLabel[])
	reg = PNML.IDRegistry()
	pntd = PnmlCore()
	PNML.add_label!(d,root(EzXML.parsexml("<test1> 1 </test1>")), pntd; reg)
	PNML.add_label!(d,root(EzXML.parsexml("<test2> 2.0 </test2>")), pntd; reg)
	PNML.add_label!(d,root(EzXML.parsexml("<test3> true </test3>")), pntd; reg)
	PNML.add_label!(d,root(EzXML.parsexml("<test4> Pnml.Term() </test4>")), pntd; reg)
end

# ╔═╡ ddacb450-31cb-4f7d-b62f-c81d0c2df858
md"""
## Step 1: Define the states and transitions of the Petri Net

Here we have 2 states, wolves and rabbits, and transitions to model predation between the two species in the system.
"""

# ╔═╡ 22a06164-6096-4a7e-a205-c1686c4b03d8
xml = """<?xml version="1.0"?>
    <pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
        <net id="net0" type="nonstandard">
        <page id="page0">
            <place id="rabbits"> <initialMarking> <text>100.0</text> </initialMarking> </place>
            <place id="wolves">  <initialMarking> <text>10.0</text> </initialMarking> </place>
            <transition id ="birth">     <rate> <text>0.3</text> </rate> </transition>
            <transition id ="death">     <rate> <text>0.7</text> </rate> </transition>
            <transition id ="predation"> <rate> <text>0.015</text> </rate> </transition>
            <arc id="a1" source="rabbits"   target="birth">     <inscription><text>1</text> </inscription> </arc>
            <arc id="a2" source="birth"     target="rabbits">   <inscription><text>2</text> </inscription> </arc>
            <arc id="a3" source="wolves"    target="predation"> <inscription><text>1</text> </inscription> </arc>
            <arc id="a4" source="rabbits"   target="predation"> <inscription><text>1</text> </inscription> </arc>
            <arc id="a5" source="predation" target="wolves">    <inscription><text>2</text> </inscription> </arc>
            <arc id="a6" source="wolves"    target="death">     <inscription><text>1</text> </inscription> </arc>
        </page>
        </net>
    </pnml>
"""

# ╔═╡ 706313c3-946c-46cd-a32d-05a4103ba9a4
md"""
## Step 2:Define the parameters and transition rates

Once a model is defined, we can define out initial parameters `u0`, a time
span `tspan`, and the transition rates of the interactions `β`.
"""

# ╔═╡ c8bf4d0a-f9a2-4a67-85f1-16ff4e53eea2
tspan = (0.0,100.0)

# ╔═╡ a3044403-467e-4002-9139-69dfefbf9990
md"""
## Step 3: Generate a solver and solve
Finally we can generate a solver and solve the simulation.
"""

# ╔═╡ 5979a601-8b07-4f33-a766-f8e3cd421891
md"""
birth rate $(@bind BR_slider Slider(0.0:0.01:1.0; show_value=true, default=0.3))
"""

# ╔═╡ edcda833-758c-4405-9459-0ba6b8c5dd3f

# ╔═╡ db25d962-4bb3-44b2-b844-ae4e4b8a9f37
md"""
death rate $(@bind DR_slider Slider(0.0:0.01:1.0; show_value=true, default=0.7))
"""

# ╔═╡ 502fbe5a-8ef5-42d6-84b1-3bf09d8c4998
md"""
predatation rate $(@bind PR_slider Slider(0.0:0.001:0.2; show_value=true, default=0.015))
"""

# ╔═╡ 55da1e63-bdef-4eff-8353-48a2a0e4bb1c
# Experimental versions on the XML.
xml2 = """<?xml version="1.0"?>
<pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
<net id="net0" type="continuous">
<page id="page0">
  <place id="rabbits"> <initialMarking>100.0</initialMarking> </place>
  <place id="wolves">  <initialMarking> <text>10.0</text> </initialMarking> </place>
  <transition id ="birth">     <rate> $BR_slider </rate> </transition>
  <transition id ="death">     <rate> $DR_slider </rate> </transition>
  <transition id ="predation"> <rate> $PR_slider </rate> </transition>
  <arc id="a1" source="rabbits"   target="birth">     <inscription><text>1</text> </inscription> </arc>
  <arc id="a2" source="birth"     target="rabbits">   <inscription><text>2</text> </inscription> </arc>
  <arc id="a3" source="wolves"    target="predation"> <inscription><text>1</text> </inscription> </arc>
  <arc id="a4" source="rabbits"   target="predation"> <inscription><text>1</text> </inscription> </arc>
  <arc id="a5" source="predation" target="wolves">    <inscription><text>2</text> </inscription> </arc>
  <arc id="a6" source="wolves"    target="death">     <inscription><text>1</text> </inscription> </arc>
</page>
</net>
</pnml>
"""

# ╔═╡ fa79b9c4-188a-4cde-9af3-b956cb0a60cd
net = PNML.SimpleNet(xml2)

# ╔═╡ a73fda08-7414-4d21-9f0e-0edf5f7cb876
# ╠═╡ show_logs = false
lotka = let
	# **Step 1:**
	@show S = PNML.place_ids(net) # [:rabbits, :wolves]
	@show Δ = PNML.transition_function(net)
	Petri.Model(S, Δ)
end;

# ╔═╡ c1dbe979-b219-452f-9f4c-0e4de84d4395
Graph(lotka)

# ╔═╡ 0c300670-010b-42c9-8f62-f681668050c3
u0 = PNML.initialMarking(net)

# ╔═╡ 2a1dc647-9aa8-4907-89f3-a15f13b45eb2
β = PNML.rates(net)

# ╔═╡ 3bc171e2-307d-4868-bd7a-e44857a1131e
prob = ODEProblem(lotka, u0, tspan, β)

# ╔═╡ 512a12f5-fd88-432b-9a56-880e2bf6879a
sol = OrdinaryDiffEq.solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8);

# ╔═╡ 91f5cd89-8c4b-4ffb-958e-25dafebbf77d
plot(sol)

# ╔═╡ Cell order:
# ╟─e6e05790-dd20-11ec-39b0-eb0589d76e44
# ╠═0143155c-a570-4a03-9cab-5704b7e49968
# ╠═a68432ac-b7eb-4985-8232-67bde7bcfd0a
# ╠═cff98b32-c6d1-4c4d-beba-ab46f9a54172
# ╠═97ed66fd-e8aa-4291-81f2-c4981c521862
# ╠═9420e502-ed69-443b-944b-9339f2e699b1
# ╟─5db3c534-07dc-4a9d-8782-41423da2ccfe
# ╠═a995864c-937e-4a0d-b99a-31fd134e23dc
# ╟─ddacb450-31cb-4f7d-b62f-c81d0c2df858
# ╟─22a06164-6096-4a7e-a205-c1686c4b03d8
# ╠═55da1e63-bdef-4eff-8353-48a2a0e4bb1c
# ╠═fa79b9c4-188a-4cde-9af3-b956cb0a60cd
# ╠═a73fda08-7414-4d21-9f0e-0edf5f7cb876
# ╠═c1dbe979-b219-452f-9f4c-0e4de84d4395
# ╟─706313c3-946c-46cd-a32d-05a4103ba9a4
# ╠═0c300670-010b-42c9-8f62-f681668050c3
# ╠═2a1dc647-9aa8-4907-89f3-a15f13b45eb2
# ╠═c8bf4d0a-f9a2-4a67-85f1-16ff4e53eea2
# ╟─a3044403-467e-4002-9139-69dfefbf9990
# ╠═3bc171e2-307d-4868-bd7a-e44857a1131e
# ╠═512a12f5-fd88-432b-9a56-880e2bf6879a
# ╠═5979a601-8b07-4f33-a766-f8e3cd421891
# ╠═edcda833-758c-4405-9459-0ba6b8c5dd3f
# ╠═db25d962-4bb3-44b2-b844-ae4e4b8a9f37
# ╠═502fbe5a-8ef5-42d6-84b1-3bf09d8c4998
# ╠═91f5cd89-8c4b-4ffb-958e-25dafebbf77d