text:
EXAMPLES_DIR = "../examples/"

function run_example(replace_strings, example_name, module_suffix="")
    example_filepath = joinpath(EXAMPLES_DIR, example_name * ".jl")
    txt = read(example_filepath, String)

    for strs in replace_strings
        txt = replace(txt, strs[1] => strs[2])
    end

    test_script_filepath = example_name * "_example_test.jl"

    open(test_script_filepath, "w") do f
        write(f, "module Test_$example_name" * "_$module_suffix\n")
        write(f, txt)
        write(f, "\nend # module")
    end

    try
        include(test_script_filepath)
    catch err
        @error sprint(showerror, err)
        rm(test_script_filepath)
        return false
    end

    rm(test_script_filepath)
    return true
end

@testset "Examples" begin
    println("Testing examples...")

    for arch in archs
        @testset "Wind and convection mixing example [$(typeof(arch))]" begin
            println("  Testing wind and convection-driving mixing example [$(typeof(arch))]")

            replace_strings = [ ("Nz = 48", "Nz = 16"),
                                ("while model.clock.time < end_time", "while model.clock.iteration < 1"),
                                ("time_step!(model, 10, wizard.Δt)", "time_step!(model, 1, wizard.Δt)"),
                              ]

            arch == GPU() && push!(replace_strings, ("architecture = CPU()", "architecture = GPU()"))

            @test_skip run_example(replace_strings, "ocean_wind_mixing_and_convection", string(typeof(arch)))
            rm("ocean_wind_mixing_and_convection.jld2", force=true)
        end
    end

    @testset "Ocean convection with plankton example" begin
        println("  Testing ocean convection with plankton example")

        replace_strings = [ ("Nz = 128", "Nz = 16"),
                           ("end_time = 1day", "end_time=1.0"),
                           ("time_step!(model, 10", "time_step!(model, 1")
                          ]

        @test_skip run_example(replace_strings, "ocean_convection_with_plankton")
    end

    @testset "Simple diffusion example" begin
        println("  Testing simple diffusion example")

        replace_strings = [ ("N = (1, 1, 128)", "N = (1, 1, 16)"),
                            ("Nt = 1000", "Nt = 2")
                          ]

        @test_skip run_example(replace_strings, "simple_diffusion")
    end

    @testset "Internal wave example" begin
        println("  Testing internal wave example")

        replace_strings = [ ("Nx = 128", "Nx = 16"),
                            ("i = 1:10", "i = 1:1"),
                            ("Nt = 200", "Nt = 2"),
                          ]

        @test_skip run_example(replace_strings, "internal_wave")
    end

    @testset "Two-dimensional turbulence example" begin
        println("  Testing two-dimensional turbulence example")

        replace_strings = [ ("N=(128, 128, 1)", "N=(16, 16, 1)"),
                            ("i = 1:10", "i = 1:1"),
                            ("Nt = 100", "Nt = 2")
                          ]

        @test_skip run_example(replace_strings, "two_dimensional_turbulence")
    end

end


text:
using Test, Colors

@testset "Colormaps" begin
    @testset "distinguishable_colors" begin
        cols = distinguishable_colors(10)
        @test isconcretetype(eltype(cols))

        mindiff = Inf
        for i = 1:10, j = i+1:10
            mindiff = min(mindiff, colordiff(cols[i], cols[j]))
        end
        @test mindiff > 8

        seed = distinguishable_colors(1)
        @test colordiff(distinguishable_colors(1, seed)[1], seed[1]) == 0.0
        @test colordiff(distinguishable_colors(1, seed; dropseed=true)[1], seed[1]) > 50

        cols_i = distinguishable_colors(4, LCHab(60, 50, 40), transform=identity,
                                        lchoices=[40, 60], cchoices=[50], hchoices=[40, 140])
        cols_p = distinguishable_colors(4, LCHab(60, 50, 40), transform=protanopic,
                                        lchoices=[40, 60], cchoices=[50], hchoices=[40, 140])
        @test cols_i[2] ≈ LCHab(40, 50, 140) atol=1f-3 # green
        @test cols_p[2] ≈ LCHab(40, 50,  40) atol=1f-3 # red
    end

    @test length(colormap("RdBu", 100)) == 100

    @test isconcretetype(eltype(colormap("Grays")))

    @test_throws ArgumentError colormap("Grays", N=10) # optional arguments, not keyword

    # not return values to check here, just checking that keywords can be used
    # Sequential
    default_blues = colormap("Blues", 10)
    @test colormap("Blues", 10; mid=0.9) == default_blues # `mid` is for diverging colormaps only (issue #300)
    @test colormap("Blues", 10; logscale=true) != default_blues
    @test colormap("Blues", 10; h=0.5) != default_blues
    @test colormap("Blues", 10; w=0.5) != default_blues
    @test colormap("Blues", 10; d=0.5) != default_blues
    @test colormap("Blues", 10; c=0.5) != default_blues
    @test colormap("Blues", 10; s=0.5) != default_blues
    @test colormap("Blues", 10; b=0.5) != default_blues
    @test colormap("Blues", 10; wcolor=colorant"white") != default_blues
    @test colormap("Blues", 10; dcolor=colorant"black") != default_blues
    @test_throws ArgumentError colormap("Blues", 10; h1=0.5)

    # Diverging
    default_rdbu = colormap("RdBu", 10)
    @test colormap("RdBu", 10; mid=0.9) != default_rdbu
    @test colormap("RdBu", 10; logscale=true) != default_rdbu
    @test colormap("RdBu", 10; h1=0.5) != default_rdbu
    @test colormap("RdBu", 10; h2=0.5) != default_rdbu
    @test colormap("RdBu", 10; w=0.5) != default_rdbu
    @test colormap("RdBu", 10; d1=0.5) != default_rdbu
    @test colormap("RdBu", 10; d2=0.5) != default_rdbu
    @test colormap("RdBu", 10; c=0.5) != default_rdbu
    @test colormap("RdBu", 10; s=0.5) != default_rdbu
    @test colormap("RdBu", 10; b=0.5) != default_rdbu
    @test colormap("RdBu", 10; wcolor=colorant"white") != default_rdbu
    @test colormap("RdBu", 10; dcolor1=colorant"black") != default_rdbu
    # `dcolor2` is disabled by the default `d2`(=0)
    @test colormap("RdBu", 10; dcolor2=colorant"black", d2=1) != default_rdbu
    @test_throws ArgumentError colormap("RdBu", 10; h=0.5)

    # The outputs of `colormap()` were slightly affected by the bug fix of
    # `MSC(h)` (issue #349).
    # The following were generated by `colormap()` in Colors.jl v0.9.6.
    blues_old   = ["F4FDFF", "B3E3F4", "65B9E7", "2978BE", "0B2857"]
    greens_old  = ["FAFFF7", "B6EEA0", "75C769", "308B40", "00391A"]
    grays_old   = ["FFFFFF", "DCDCDC", "A9A9A9", "626262", "000000"]
    oranges_old = ["FFFBF6", "FFD6B4", "FF9D5F", "E2500D", "732108"]
    purples_old = ["FBFBFB", "DBD7F6", "AFA7E6", "7666BF", "3C0468"]
    reds_old    = ["FFF1EE", "FFC4B9", "FF8576", "E72823", "6D0B0C"]
    rdbu_old    = ["610102", "FF8D7B", "F9F8F9", "76B4E8", "092C58"]

    to_rgb(s) = parse(RGB, "#"*s)
    max_colordiff(a1, a2) = reduce(max, colordiff.(a1, a2))
    @test max_colordiff(colormap("Blues", 5), to_rgb.(blues_old)) < 1
    @test max_colordiff(colormap("Greens", 5), to_rgb.(greens_old)) < 1
    @test max_colordiff(colormap("Grays", 5), to_rgb.(grays_old)) < 1
    @test max_colordiff(colormap("Oranges", 5), to_rgb.(oranges_old)) < 1
    @test max_colordiff(colormap("Purples", 5), to_rgb.(purples_old)) < 1
    @test max_colordiff(colormap("Reds", 5), to_rgb.(reds_old)) < 1
    @test max_colordiff(colormap("RdBu", 5), to_rgb.(rdbu_old)) < 1

    # TODO: add more tests
end

ext: jl


text:
# ---
# title: 199. Binary Tree Right Side View
# id: problem199
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/binary-tree-right-side-view/description/>
# hidden: true
# ---
#
# Given a binary tree, imagine yourself standing on the _right_ side of it,
# return the values of the nodes you can see ordered from top to bottom.
#
# **Example:**
#
#
#
#     Input:  [1,2,3,null,5,null,4]
#     Output:  [1, 3, 4]
#     Explanation:
#        1            <---
#      /   \
#     2     3         <---
#      \     \
#       5     4       <---
#
#
#
## @lc code=start
using LeetCode

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


text:
#= This script includes all modules and functions, sets up the grid parameters
etc, and runs the model =#

using Dates
using NetCDF
using FileIO

if VERSION == v"0.7.0"
    using Printf
    #using Statistics
end

#using MPI
#using SigmoidNumbers

# Finite16nonu
#include("/home/kloewer/julia/FiniteFloats.jl/src/FiniteFloats.jl")

# PARAMETERS, GRID, CONSTANTS and DOMAIN DECOMPOSITION
include("parameters.jl")
include("src/grid.jl")
include("src/constants.jl")
#include("src/domain_decomposition.jl")

# OPERATORS and everything that is needed for the RHS
include("src/gradients.jl")
include("src/interpolations.jl")
include("src/PV_adv.jl")
include("src/bottom_friction.jl")
include("src/diffusion.jl")
include("src/tracer_adv.jl")
include("src/coriolis.jl")
include("src/forcing.jl")
include("src/bottom_topography.jl")
include("src/rhs.jl")
include("src/time_integration.jl")
include("src/ghost_points.jl")
include("src/initial_conditions.jl")
include("src/preallocate.jl")

# OUTPUT AND FEEDBACK
include("src/feedback.jl")
include("src/output.jl")
global run_id,runpath
run_id,runpath = get_run_id_path()

# INITIALISE & RUN
u,v,η,sst = initial_conditions()
u,v,η,sst = time_integration(u,v,η,sst)

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
# Magnetic dipole field.

"Analytic electric field function for testing."
function getE_dipole(xu)
   [0.0, 0.0, 0.0]
end

"Analytic magnetic field function for testing. Return in SI unit."
function getB_dipole(xu)
   BMoment = BMoment_Earth
   dipole(xu[1:3], BMoment)
end

"Calculates the magnetic field from a dipole with magnetic moment `M` at `r`."
function dipole(rIn, M)
   x, y, z = rIn
   r = sqrt(x^2 + y^2 + z^2)
   Coef = μ₀/(4*π*r^5)

   B = [3*x^2-r^2 3*x*y     3*x*z;
        3*y*x     3*y^2-r^2 3*y*z;
        3*z*x     3*z*y     3*z^2-r^2
       ] * M * Coef
end

"""
    fieldline(L, ϕ, nP)

Creates points on one field line of the magnetic field from a dipole. In a centered dipole
magnetic field model, the path along a given L shell can be described as r = L*cos²λ,
where r is the radial distance (in planetary radii) to a point on the line,
λ is its co-latitude, and L is the L-shell of interest.
"""
function fieldline(ϕ::Float64, L::Float64=2.5, nP::Int=100)

   xyz = [ sph2cart(L*sin(θ)^2,ϕ,θ) for θ in range(-π,stop=π,length=nP) ]
   x = Vector{Float64}(undef,length(xyz))
   y = Vector{Float64}(undef,length(xyz))
   z = Vector{Float64}(undef,length(xyz))

   for (i, pos) in enumerate(xyz)
      x[i],y[i],z[i] = [pos...]
   end

   (x,y,z)
end


text:
# Autogenerated wrapper script for CeresCApi_jll for x86_64-linux-gnu-cxx11
export libceres_c_api

using CeresSolver_jll
JLLWrappers.@generate_wrapper_header("CeresCApi")
JLLWrappers.@declare_library_product(libceres_c_api, "libceres_c_api.so")
function __init__()
    JLLWrappers.@generate_init_header(CeresSolver_jll)
    JLLWrappers.@init_library_product(
        libceres_c_api,
        "lib/libceres_c_api.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using QuasiHamiltonianRTN
using QuasiHamiltonianRTN.Lattice1D
using Combinatorics
using SparseArrays
using LinearAlgebra

include("../src/SpecialUnitary.jl")
include("../src/Utilities.jl")
using Test

@testset "Utilities" begin
    @testset "random_density_matrix" for i = 1:10
        ρ = random_density_matrix(4)
        @test tr(ρ) ≈ 1
        @test ρ ≈ ρ'
        @test real(tr(ρ^2)) <= 1.
        @test all(eigvals(ρ) .>= 0)
    end
end

@testset "SU(N) generators" begin
    @testset "N = $N" for N = 20:21
        let λ = sun_generators(N)
            @testset "Traceless" begin
                @test mapreduce(A -> isapprox(tr(A), 0., atol=eps(Float64) * N), &, λ)
            end
            # Hermitian
            @testset "Hermitian" begin
                @test mapreduce(A -> A ≈ A', &, λ)
            end
            # Orthogonal
            @testset "Orthogonal" begin
                @test [tr(λ[i] * λ[j]) for i = 1:length(λ), j = 1:length(λ)] ≈ 2*I
             end
        end
    end
end

@testset "SU(N) structure constants" begin
    @testset "SU(2) gives the Levi-Civita tensor" begin
        @test reduce(&, [structure_constants(sun_generators(2))[i][j,k] == levicivita([i,j,k]) for i=1:3,j=1:3,k=1:3])
    end

    @testset "[λ[i], λ[j]] = i ∑ f[i][j,k] λ[k] for SU($n)" for n = 7:9
        let λ = sun_generators(n), f = structure_constants(λ)
            @test reduce(&,[commutator(λ[i], λ[j]) ≈ 2*im*sum([f[i][j,k]*λ[k] for k = 1:n^2-1])  for i = 1 : n^2-1, j = 1:n^2-1])
        end
    end
end

@testset "Bloch vector" begin
    @testset "Back and forth" begin
        for n = 2:10
            ρ = random_density_matrix(n)
            @test Matrix(density_operator(bloch_vector(ρ))) ≈ ρ
        end
    end
end

@testset "Unitary evolution check" begin
    let n = 4
        ρ = random_density_matrix(n)
        ρ1 = unitary_op(Lattice1D.H0(n), 1.)*ρ*unitary_op(Lattice1D.H0(n), 1.)'
        ρ2 = Matrix(density_operator(exp(- Matrix(igen(Lattice1D.H0(n)))) * bloch_vector(ρ)))
        @test ρ1 ≈ ρ2
    end
end

@testset "Noiseless evolution check" begin
    let tv = range(0, stop=5, length=10), n = 5, ρ0 = localized_state(n, 3)
        rtU = [bloch_vector(unitary_op(Lattice1D.H0(n), t) * ρ0 * unitary_op(Lattice1D.H0(n), t)') for t in tv]
        rtE = evolution(Lattice1D.Hamiltonian(n, .0), bloch_vector(ρ0), tv, γ=0.)
        @test hcat(rtU...) ≈ rtE

        Hq = quasiHamiltonian(Lattice1D.Hamiltonian(n, .0,), 0.)
        rtEHq = evolution(Hq, bloch_vector(ρ0), tv)

        @test rtEHq ≈ rtE

        rtEHqEK = evolution(Hq, bloch_vector(ρ0), tv; expmpkg=:Expokit)
        @test rtEHq ≈ rtEHqEK
    end
end

@testset "Noisy evolution" begin
    let t = 1., longt = 80, n = 5, ρ0 = localized_state(n, 3)
        H = Lattice1D.Hamiltonian(n, 1.)
        Hq = quasiHamiltonian(H, 1.)

        # Check that the two expm give similar results
        rtEHq = evolution(Hq, bloch_vector(ρ0), t)
        rtEHqEK = evolution(Hq, bloch_vector(ρ0), t; expmpkg=:Expokit)

        @test rtEHq ≈ rtEHqEK

        # Check that the stationary state is the maximally mixed one
        # (i.e. Bloch vector is zero)

        rtEHq = evolution(Hq, bloch_vector(ρ0), longt)
        rtEHqEK = evolution(Hq, bloch_vector(ρ0), longt; expmpkg=:Expokit)
        @test rtEHq ≈ zero(rtEHq) atol=1e-16
    end
end


text:
using VibronicToolkit: jackknife

using Random: MersenneTwister
using Statistics: cov, mean, std, var

# This uses the formulas from Peter Young's notes
# (https://arxiv.org/abs/1210.3781, https://arxiv.org/pdf/1210.3781v3.pdf).
function propagate_errors_nonlinear1(f, dfx, dfxx, xs)
    mu_x = mean(xs)
    biased_mean = f(mu_x)
    pre_bias = dfxx(mu_x)*var(xs)
    unbiased_mean = biased_mean - 0.5*pre_bias/length(xs)
    pre_err = dfx(mu_x)^2*var(xs)
    err = sqrt(pre_err/length(xs))
    unbiased_mean, err
end

function propagate_errors_nonlinear2(f, dfx, dfy, dfxx, dfxy, dfyy, xs, ys)
    mu_x, mu_y = mean(xs), mean(ys)
    biased_mean = f(mu_x,mu_y)
    pre_bias = dfxx(mu_x,mu_y)*var(xs) + 2dfxy(mu_x,mu_y)*cov(xs,ys) + dfyy(mu_x,mu_y)*var(ys)
    unbiased_mean = biased_mean - 0.5*pre_bias/length(xs)
    pre_err = dfx(mu_x,mu_y)^2*var(xs) + 2dfx(mu_x,mu_y)*dfy(mu_x,mu_y)*cov(xs,ys) + dfy(mu_x,mu_y)^2*var(ys)
    err = sqrt(pre_err/length(xs))
    unbiased_mean, err
end

compare_stats((mu1, err1), (mu2, err2)) = isapprox(mu1, mu2; rtol=1e-11) && isapprox(err1, err2; rtol=1e-6)

let rng = MersenneTwister(0x1234),
    xs = randn(rng, 10000),
    ys = randn(rng, 10000)

    # Linear function of a single variable.
    l1a(x) = 2x.-3.7
    @test compare_stats(jackknife(l1a, xs),
                        (l1a(mean(xs)),
                         2std(xs)/sqrt(length(xs))))
    l1b(y) = -1.8y.+9
    @test compare_stats(jackknife(l1b, ys),
                        (l1b(mean(ys)),
                         1.8std(ys)/sqrt(length(ys))))

    # Linear function of two variables.
    l2a(x, y) = l1a(x) + l1b(y)
    @test compare_stats(jackknife(l2a, xs, ys),
                        (l2a(mean(xs), mean(ys)),
                         sqrt((2std(xs))^2+(1.8std(ys))^2-2*2*1.8*cov(xs,ys))/sqrt(length(xs))))

    # Nonlinear function of one variable.
    n1a(y) = 1 ./ l1b(y)
    dn1ay(y) = 1.8 ./ l1b(y)^2
    dn1ayy(y) = 2*1.8^2 ./ l1b(y)^3
    @test compare_stats(jackknife(n1a, ys),
                        (propagate_errors_nonlinear1(n1a, dn1ay, dn1ayy, ys)))

    # Nonlinear function of two variables.
    n2a(x, y) = l1a(x) ./ l1b(y)
    dn2ax(x, y) = 2 ./ l1b(y)
    dn2ay(x, y) = 1.8*l1a(x) ./ l1b(y)^2
    dn2axx(x, y) = 0.0
    dn2axy(x, y) = 2*1.8 ./ l1b(y)^2
    dn2ayy(x, y) = 2*1.8^2*l1a(x) ./ l1b(y)^3
    @test compare_stats(jackknife(n2a, xs, ys),
                        propagate_errors_nonlinear2(n2a, dn2ax, dn2ay, dn2axx, dn2axy, dn2ayy, xs, ys))

    @test_throws DomainError jackknife(n2a)
    @test_throws DomainError jackknife(n2a, [1.0], [2.0, 3.0])
end


text:
using EngEconomics, Roots, Plots, ColorSchemes
plotly()

# Given
# - Finedetail = A, Simplicity = B
strA = "Finedetail"
strB = "Simplicity"
lifeA = 7
lifeB = 10
PA = 200000 # Initial Cost of A - [$]
PB = 350000 # Initial Cost of B - [$]
MA = 10000 # Maintenance Cost of A - [$/yr]
UMA = 0.05 # Maintenance Cost of A - [$/units]
MB = 20000 # Maintenance Cost of B - [$/yr]
UMB = 0.01 # Maintenance Cost of B - [$/units]
LA = 1.25 # Labor Cost of A - [$/unit]
LB = 0.50 # Labor Cost of B - [$/unit]
OA = 6500 # "Other" Costs of A - [$/yr]
UOA = 0.95 # "Other" Costs of A - [$/unit]
OB = 15500 # "Other" Costs of B - [$/yr]
UOB = 0.55 # "Other" Costs of B - [$/unit]
FA = 5000 # Salvage Value of A - [$]
FB = 20000 # Salvage Value of B - [$]
MARR = 0.15

# Find
# (a.) Who is the preferred supplier if sales are 30,000 units per year? Use an
#      annual worth comparison.
unitRate = 30000 # [units/yr]
EUAWA = -PA * capitalRecoveryFactor(MARR, lifeA) + -MA + -UMA * unitRate + -LA * unitRate + -OA + -UOA * unitRate + FA * sinkingFundFactor(MARR, lifeA)
EUAWB = -PB * capitalRecoveryFactor(MARR, lifeB) + -MB + -UMB * unitRate + -LB * unitRate + -OB + -UOB * unitRate + FB * sinkingFundFactor(MARR, lifeB)

if EUAWA > EUAWB
	println(strA, " is the preferred supplier if unit sales are ", unitRate, " because its EUAW exceeds that of ", strB)
else
	println(strB, " is the preferred supplier if unit sales are ", unitRate, " because its EUAW exceeds that of ", strA)
end

# (b.) Who is the preferred supplier if sales are 200,000 units per year. Use
#      annual worth comparison
unitRate = 200000 # [units/yr]
EUAWA = -PA * capitalRecoveryFactor(MARR, lifeA) + -MA + -UMA * unitRate + -LA * unitRate + -OA + -UOA * unitRate + FA * sinkingFundFactor(MARR, lifeA)
EUAWB = -PB * capitalRecoveryFactor(MARR, lifeB) + -MB + -UMB * unitRate + -LB * unitRate + -OB + -UOB * unitRate + FB * sinkingFundFactor(MARR, lifeB)

if EUAWA > EUAWB
	println(strA, " is the preferred supplier if unit sales are ", unitRate, " because its EUAW exceeds that of ", strB)
else
	println(strB, " is the preferred supplier if unit sales are ", unitRate, " because its EUAW exceeds that of ", strA)
end

# (c.) How sensitive is the choice of supplier to sales level? Experiment with
#      sales levels between 30,000 and 200,000 units per year. At what sales
#      level will the costs of the two melters be equal?
unitRateVector = collect(30000:1000:200000)
EUAWA = -PA * capitalRecoveryFactor(MARR, lifeA) .+ -MA .+ -UMA * unitRateVector .+ -LA * unitRateVector .+ -OA .+ -UOA * unitRateVector .+ FA * sinkingFundFactor(MARR, lifeA)
EUAWB = -PB * capitalRecoveryFactor(MARR, lifeB) .+ -MB .+ -UMB * unitRateVector .+ -LB * unitRateVector .+ -OB .+ -UOB * unitRateVector .+ FB * sinkingFundFactor(MARR, lifeB)

plot(unitRateVector, EUAWA, xlabel="Units", ylabel="Annual Worth (\$)", label=strA,
	title="Sensitivity of Alternatives to Unit Sales", legend=:outerbottomright)
plot!(unitRateVector, EUAWB, label=strB)

EUAWA_func(u) = -PA * capitalRecoveryFactor(MARR, lifeA) .+ -MA .+ -UMA * u .+ -LA * u .+ -OA .+ -UOA * u .+ FA * sinkingFundFactor(MARR, lifeA)
EUAWB_func(u) = -PB * capitalRecoveryFactor(MARR, lifeB) .+ -MB .+ -UMB * u .+ -LB * u .+ -OB .+ -UOB * u .+ FB * sinkingFundFactor(MARR, lifeB)
f(u) = EUAWA_func(u) - EUAWB_func(u)
unitBreakEven = find_zero(f, 34000)

println("The annual worths of both alternatives are sensitive to the unit sales rate.")
println("The breakeven unit sales rate is: ", unitBreakEven)
println("From the plot, we understand that if sales are below the breakeven unit sales rate, then Finedetail is preferred. However, if the unit sales rate exceeds the breakeven unit sales rate, then simplicity is preferred")


text:
module ComplexOperations

using StaticArrays

export conjugate

export complex_multiply
export complex_multiply_SMatrix

export complex_divide
export complex_divide_SMatrix

export complex_vector_dot_product
export complex_vector_dot_product_SMatrix

export complex_vector_cross_product
export complex_vector_cross_product_SMatrix

export complex_matrix_inversion
export complex_matrix_multiplication

# Conjugate -----------------------------------------------------------------------------------
"""
    Calculate conjugate of imaginary number
"""
function conjugate(a::R, b::R) where R <: Real
    return hcat(a, -b)
end

function conjugate(a::AbstractVecOrMat{R}, b::AbstractVecOrMat{R}) where R <: Real
    return hcat(a, -b)
end

# Scalar multiplication -----------------------------------------------------------------------------------
"""
    Multiplication of two complex numbers Z1=a+im*b and Z2=c+im*d, defined by a,b,c and d.
a, b, c, and d are all real numbers.
Multiplication algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_multiply(a::R, b::R, c::R, d::R) where R <: Real
    return @fastmath @inbounds hcat(a * c - b * d, b * c + a * d)
end

"""
    Multiplication of two vectors Z1=a+im*b and Z2=c+im*d.
a, b, c, and d are all arrays of real numbers, all arrays have the same size.
Multiplication algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_multiply(a::AbstractVecOrMat{R}, b::AbstractVecOrMat{R}, c::AbstractVecOrMat{R}, d::AbstractVecOrMat{R}) where R <: Real
    return @fastmath @inbounds hcat(a .* c - b .* d, b .* c + a .* d)
end

"""
    Multiplication of two vectors Z1 and Z2.
We assume the first and second columns of each vector are the real and imaginary parts, respectively
Each of Z1 and Z2 has two columns, and arbitrary number of rows.
Multiplication algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_multiply(Z1::AbstractVecOrMat{R}, Z2::AbstractVecOrMat{R}) where R <: Real
    return vcat(complex_multiply.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Scalar multiplication, returns SMatrix -----------------------------------------------------------------------------------
"""
    Same as `complex_multiply`, but return SMatrix
"""
function complex_multiply_SMatrix(a::R, b::R, c::R, d::R) where R <: Real
    return @fastmath @inbounds SMatrix{1,2}(a * c - b * d, b * c + a * d)
end

# TODO: remove this one, it is confusing
"""
    Same as `complex_multiply`, but return SMatrix
"""
function complex_multiply_SMatrix(a::AbstractVecOrMat{R}, b::AbstractVecOrMat{R}, c::AbstractVecOrMat{R}, d::AbstractVecOrMat{R}) where R <: Real
    return @fastmath @inbounds vcat(complex_multiply_SMatrix.(a, b, c, d)...)
end

"""
    Same as `complex_multiply`, but return SMatrix
"""
function complex_multiply_SMatrix(Z1::AbstractVecOrMat{R}, Z2::AbstractVecOrMat{R}) where R <: Real
    return vcat(complex_multiply_SMatrix.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Scalar division -----------------------------------------------------------------------------------
"""
    Division of two vectors Z1=a+im*b and Z2=c+im*d.
a, b, c, and d are all real numbers
Division algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_divide(a::R, b::R, c::R, d::R) where R <: Real
    return @fastmath @inbounds hcat((a * c + b * d) / (c^2 + d^2), (b * c - a * d) / (c^2 + d^2))
end

"""
    Division of two vectors Z1=a+im*b and Z2=c+im*d.
a, b, c, and d are all arrays of real numbers, all arrays have the same size.
Division algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_divide(a::AbstractVecOrMat{R}, b::AbstractVecOrMat{R}, c::AbstractVecOrMat{R}, d::AbstractVecOrMat{R}) where R <: Real
    return @fastmath @inbounds hcat((a .* c + b .* d) / (c.^2 + d.^2), (b .* c - a .* d) / (c.^2 + d.^2))
end

"""
    Division of two vectors Z1 and Z2.
We assume the first and second columns of each vector are the real and imaginary parts, respectively
Each of Z1 and Z2 has two columns, and arbitrary number of rows.
Division algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_divide(Z1::AbstractVecOrMat{R}, Z2::AbstractVecOrMat{R}) where R <: Real
    return vcat(complex_divide.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Scalar division, returns SMatrix -----------------------------------------------------------------------------------
"""
    Same as `complex_divide`, but returns SMatrix
"""
function complex_divide_SMatrix(a::R, b::R, c::R, d::R) where R <: Real
    return @fastmath @inbounds SMatrix{1,2}((a * c + b * d) / (c^2 + d^2), (b * c - a * d) / (c^2 + d^2))
end

# TODO: remove this one, it is confusing
"""
    Same as `complex_divide`, but returns SMatrix
"""
function complex_divide_SMatrix(a::AbstractVecOrMat{R}, b::AbstractVecOrMat{R}, c::AbstractVecOrMat{R}, d::AbstractVecOrMat{R}) where R <: Real
    return @fastmath @inbounds vcat(complex_divide_SMatrix.(a, b, c, d)...)
end
"""
    Same as `complex_divide`, but returns SMatrix
"""
function complex_divide_SMatrix(Z1::AbstractVecOrMat{R}, Z2::AbstractVecOrMat{R}) where R <: Real
    return vcat(complex_divide_SMatrix.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Vector dot product ------------------------------------------------------------------------------------------
"""
    Claculate dot product of two vectors
Inputs are real and imaginary parts of the two vectors, each is represented as 3-element array.
"""
function complex_vector_dot_product(A_r::AbstractVector{R}, A_i::AbstractVector{R}, B_r::AbstractVector{R}, B_i::AbstractVector{R}) where R <: Real
    return (
    complex_multiply(A_r[1], A_i[1], B_r[1], B_i[1]) +
    complex_multiply(A_r[2], A_i[2], B_r[2], B_i[2]) +
    complex_multiply(A_r[3], A_i[3], B_r[3], B_i[3])
)
end

"""
    Claculate dot product of two vectors. A and B represents complex vectors.
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_dot_product(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent real and complex vectors, respectively.
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product(A::AbstractVector{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_dot_product(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent complex and real vectors, respectively.
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product(A::AbstractMatrix{R}, B::AbstractVector{R}) where R <: Real
    return complex_vector_dot_product(A[:,1], A[:,2], B, zero(B))
end

# Vector dot product, returns SMatrix ------------------------------------------------------------------------------------------
"""
    Claculate dot product of two vectors, and returns SMatrix.
Inputs are real and imaginary parts of the two vectors, each is represented as 3-element array.
"""
function complex_vector_dot_product_SMatrix(A_r::AbstractVector{R}, A_i::AbstractVector{R}, B_r::AbstractVector{R}, B_i::AbstractVector{R}) where R <: Real
    return (
    complex_multiply_SMatrix(A_r[1], A_i[1], B_r[1], B_i[1]) +
    complex_multiply_SMatrix(A_r[2], A_i[2], B_r[2], B_i[2]) +
    complex_multiply_SMatrix(A_r[3], A_i[3], B_r[3], B_i[3])
)
end

"""
    Claculate dot product of two vectors. A and B represents complex vectors. Returns SMatrix
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product_SMatrix(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_dot_product_SMatrix(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent real and complex vectors, respectively. Returns SMatrix
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product_SMatrix(A::AbstractVector{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_dot_product_SMatrix(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent complex and real vectors, respectively. Returns SMatrix
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product_SMatrix(A::AbstractMatrix{R}, B::AbstractVector{R}) where R <: Real
    return complex_vector_dot_product_SMatrix(A[:,1], A[:,2], B, zero(B))
end

# Vector cross product ------------------------------------------------------------------------------------------
"""
    Claculate cross product of two vectors, inputs are real and imaginary parts of the two vectors
Each is represented as 3-element array.
"""
function complex_vector_cross_product(A_r::AbstractVector{R}, A_i::AbstractVector{R}, B_r::AbstractVector{R}, B_i::AbstractVector{R}) where R <: Real
    return vcat(
    complex_multiply(A_r[2], A_i[2], B_r[3], B_i[3]) - complex_multiply(A_r[3], A_i[3], B_r[2], B_i[2]),
    complex_multiply(A_r[3], A_i[3], B_r[1], B_i[1]) - complex_multiply(A_r[1], A_i[1], B_r[3], B_i[3]),
    complex_multiply(A_r[1], A_i[1], B_r[2], B_i[2]) - complex_multiply(A_r[2], A_i[2], B_r[1], B_i[1]),
)
end

"""
    Claculate cross product of two vectors. A and B represents complex vectors.
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_cross_product(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent real and complex vectors, respectively.
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product(A::AbstractVector{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_cross_product(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent complex and real vectors, respectively.
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product(A::AbstractMatrix{R}, B::AbstractVector{R}) where R <: Real
    return complex_vector_cross_product(A[:,1], A[:,2], B, zero(B))
end

# Vector cross product, returns SMatrix ------------------------------------------------------------------------------------------
"""
    Claculate cross product of two vectors, inputs are real and imaginary parts of the two vectors. Returns SMatrix
Each is represented as 3-element array.
"""
function complex_vector_cross_product_SMatrix(A_r::AbstractVector{R}, A_i::AbstractVector{R}, B_r::AbstractVector{R}, B_i::AbstractVector{R}) where R <: Real
    return vcat(
    complex_multiply_SMatrix(A_r[2], A_i[2], B_r[3], B_i[3]) - complex_multiply_SMatrix(A_r[3], A_i[3], B_r[2], B_i[2]),
    complex_multiply_SMatrix(A_r[3], A_i[3], B_r[1], B_i[1]) - complex_multiply_SMatrix(A_r[1], A_i[1], B_r[3], B_i[3]),
    complex_multiply_SMatrix(A_r[1], A_i[1], B_r[2], B_i[2]) - complex_multiply_SMatrix(A_r[2], A_i[2], B_r[1], B_i[1]),
)
end

"""
    Claculate cross product of two vectors. A and B represents complex vectors. Returns SMatrix
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product_SMatrix(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_cross_product_SMatrix(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent real and complex vectors, respectively. Returns SMatrix
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product_SMatrix(A::AbstractVector{R}, B::AbstractMatrix{R}) where R <: Real
    return complex_vector_cross_product_SMatrix(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent complex and real vectors, respectively. Returns SMatrix
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product_SMatrix(A::AbstractMatrix{R}, B::AbstractVector{R}) where R <: Real
    return complex_vector_cross_product_SMatrix(A[:,1], A[:,2], B, zero(B))
end

# Complex matrix inverse ------------------------------------------------------------------------------------------
"""
    Calculate the inverse of a complex matrix A+iB, where A and B are real matrices

There are two possible versions:
Version #1: This one preallocate. Faster and use more memory
for 1M evaluations: 1.771474 seconds (14.18 M allocations: 4.510 GiB, 4.52% gc time, 6.40% compilation time)
`
function complex_matrix_inversion(A, B)
    inv_A_B = inv(A) * B
    C = inv(A + B * inv_A_B)
    D = -inv_A_B * C
    return C, D
end
`

Version #2: No preallocation. Slower and use less memory
for 1M evaluations: 2.330454 seconds (19.00 M allocations: 6.512 GiB, 2.71% gc time)
`
function complex_matrix_inversion(A, B)
    C = inv(A + B * inv(A) * B)
    D = -1 .* inv(A) * B * C
    return C, D
end
`
"""
function complex_matrix_inversion(A::AbstractMatrix{R}, B::AbstractMatrix{R}) where R <: Real
    @fastmath @inbounds inv_A_B = inv(A) * B
    @fastmath @inbounds C = inv(A + B * inv_A_B)
    @fastmath @inbounds D = -inv_A_B * C
    return C, D
end

# Complex matrix multiplication ------------------------------------------------------------------------------------------
"""
    Multiply two matrices A+iB and C+iD, where A,B,C,D are real matrices
"""
function complex_matrix_multiplication(A::AbstractMatrix{R}, B::AbstractMatrix{R}, C::AbstractMatrix{R}, D::AbstractMatrix{R}) where R <: Real
    return @fastmath @inbounds (A * C - B * D, A * D + B * C)
end

##########################################################################################
##########################################################################################
##########################################################################################
##########################################################################################
# TODO: find a better way to accomodate different types.
# in what follows, I am adding more methods for all
# Scalar multiplication -----------------------------------------------------------------------------------
"""
    Multiplication of two complex numbers Z1=a+im*b and Z2=c+im*d, defined by a,b,c and d.
a, b, c, and d are all real numbers.
Multiplication algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_multiply(a::Real, b::Real, c::Real, d::Real)
    return @fastmath @inbounds hcat(a * c - b * d, b * c + a * d)
end

"""
    Multiplication of two vectors Z1=a+im*b and Z2=c+im*d.
a, b, c, and d are all arrays of real numbers, all arrays have the same size.
Multiplication algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_multiply(a::AbstractVecOrMat{Real}, b::AbstractVecOrMat{Real}, c::AbstractVecOrMat{Real}, d::AbstractVecOrMat{Real})
    return @fastmath @inbounds hcat(a .* c - b .* d, b .* c + a .* d)
end

"""
    Multiplication of two vectors Z1 and Z2.
We assume the first and second columns of each vector are the real and imaginary parts, respectively
Each of Z1 and Z2 has two columns, and arbitrary number of rows.
Multiplication algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_multiply(Z1::AbstractVecOrMat{Real}, Z2::AbstractVecOrMat{Real})
    return vcat(complex_multiply.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Scalar multiplication, returns SMatrix -----------------------------------------------------------------------------------
"""
    Same as `complex_multiply`, but return SMatrix
"""
function complex_multiply_SMatrix(a::Real, b::Real, c::Real, d::Real)
    return @fastmath @inbounds SMatrix{1,2}(a * c - b * d, b * c + a * d)
end

# TODO: remove this one, it is confusing
"""
    Same as `complex_multiply`, but return SMatrix
"""
function complex_multiply_SMatrix(a::AbstractVecOrMat{Real}, b::AbstractVecOrMat{Real}, c::AbstractVecOrMat{Real}, d::AbstractVecOrMat{Real})
    return @fastmath @inbounds vcat(complex_multiply_SMatrix.(a, b, c, d)...)
end

"""
    Same as `complex_multiply`, but return SMatrix
"""
function complex_multiply_SMatrix(Z1::AbstractVecOrMat{Real}, Z2::AbstractVecOrMat{Real})
    return vcat(complex_multiply_SMatrix.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Scalar division -----------------------------------------------------------------------------------
"""
    Division of two vectors Z1=a+im*b and Z2=c+im*d.
a, b, c, and d are all real numbers
Division algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_divide(a::Real, b::Real, c::Real, d::Real)
    return @fastmath @inbounds hcat((a * c + b * d) / (c^2 + d^2), (b * c - a * d) / (c^2 + d^2))
end

"""
    Division of two vectors Z1=a+im*b and Z2=c+im*d.
a, b, c, and d are all arrays of real numbers, all arrays have the same size.
Division algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_divide(a::AbstractVecOrMat{Real}, b::AbstractVecOrMat{Real}, c::AbstractVecOrMat{Real}, d::AbstractVecOrMat{Real})
    return @fastmath @inbounds hcat((a .* c + b .* d) / (c.^2 + d.^2), (b .* c - a .* d) / (c.^2 + d.^2))
end

"""
    Division of two vectors Z1 and Z2.
We assume the first and second columns of each vector are the real and imaginary parts, respectively
Each of Z1 and Z2 has two columns, and arbitrary number of rows.
Division algorithm for complex numbers can be found here: https://en.wikipedia.org/wiki/Multiplication_algorithm
"""
function complex_divide(Z1::AbstractVecOrMat{Real}, Z2::AbstractVecOrMat{Real})
    return vcat(complex_divide.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Scalar division, returns SMatrix -----------------------------------------------------------------------------------
"""
    Same as `complex_divide`, but returns SMatrix
"""
function complex_divide_SMatrix(a::Real, b::Real, c::Real, d::Real)
    return @fastmath @inbounds SMatrix{1,2}((a * c + b * d) / (c^2 + d^2), (b * c - a * d) / (c^2 + d^2))
end

# TODO: remove this one, it is confusing
"""
    Same as `complex_divide`, but returns SMatrix
"""
function complex_divide_SMatrix(a::AbstractVecOrMat{Real}, b::AbstractVecOrMat{Real}, c::AbstractVecOrMat{Real}, d::AbstractVecOrMat{Real})
    return @fastmath @inbounds vcat(complex_divide_SMatrix.(a, b, c, d)...)
end
"""
    Same as `complex_divide`, but returns SMatrix
"""
function complex_divide_SMatrix(Z1::AbstractVecOrMat{Real}, Z2::AbstractVecOrMat{Real})
    return vcat(complex_divide_SMatrix.(Z1[:,1], Z1[:,2], Z2[:,1], Z2[:,2])...)
end

# Vector dot product ------------------------------------------------------------------------------------------
"""
    Claculate dot product of two vectors
Inputs are real and imaginary parts of the two vectors, each is represented as 3-element array.
"""
function complex_vector_dot_product(A_r::AbstractVector{Real}, A_i::AbstractVector{Real}, B_r::AbstractVector{Real}, B_i::AbstractVector{Real})
    return (
    complex_multiply(A_r[1], A_i[1], B_r[1], B_i[1]) +
    complex_multiply(A_r[2], A_i[2], B_r[2], B_i[2]) +
    complex_multiply(A_r[3], A_i[3], B_r[3], B_i[3])
)
end

"""
    Claculate dot product of two vectors. A and B represents complex vectors.
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product(A::AbstractMatrix{Real}, B::AbstractMatrix{Real})
    return complex_vector_dot_product(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent real and complex vectors, respectively.
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product(A::AbstractVector{Real}, B::AbstractMatrix{Real})
    return complex_vector_dot_product(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent complex and real vectors, respectively.
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product(A::AbstractMatrix{Real}, B::AbstractVector{Real})
    return complex_vector_dot_product(A[:,1], A[:,2], B, zero(B))
end

# Vector dot product, returns SMatrix ------------------------------------------------------------------------------------------
"""
    Claculate dot product of two vectors, and returns SMatrix.
Inputs are real and imaginary parts of the two vectors, each is represented as 3-element array.
"""
function complex_vector_dot_product_SMatrix(A_r::AbstractVector{Real}, A_i::AbstractVector{Real}, B_r::AbstractVector{Real}, B_i::AbstractVector{Real})
    return (
    complex_multiply_SMatrix(A_r[1], A_i[1], B_r[1], B_i[1]) +
    complex_multiply_SMatrix(A_r[2], A_i[2], B_r[2], B_i[2]) +
    complex_multiply_SMatrix(A_r[3], A_i[3], B_r[3], B_i[3])
)
end

"""
    Claculate dot product of two vectors. A and B represents complex vectors. Returns SMatrix
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product_SMatrix(A::AbstractMatrix{Real}, B::AbstractMatrix{Real})
    return complex_vector_dot_product_SMatrix(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent real and complex vectors, respectively. Returns SMatrix
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product_SMatrix(A::AbstractVector{Real}, B::AbstractMatrix{Real})
    return complex_vector_dot_product_SMatrix(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate dot product of two vectors. A and B represent complex and real vectors, respectively. Returns SMatrix
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_dot_product_SMatrix(A::AbstractMatrix{Real}, B::AbstractVector{Real})
    return complex_vector_dot_product_SMatrix(A[:,1], A[:,2], B, zero(B))
end

# Vector cross product ------------------------------------------------------------------------------------------
"""
    Claculate cross product of two vectors, inputs are real and imaginary parts of the two vectors
Each is represented as 3-element array.
"""
function complex_vector_cross_product(A_r::AbstractVector{Real}, A_i::AbstractVector{Real}, B_r::AbstractVector{Real}, B_i::AbstractVector{Real})
    return vcat(
    complex_multiply(A_r[2], A_i[2], B_r[3], B_i[3]) - complex_multiply(A_r[3], A_i[3], B_r[2], B_i[2]),
    complex_multiply(A_r[3], A_i[3], B_r[1], B_i[1]) - complex_multiply(A_r[1], A_i[1], B_r[3], B_i[3]),
    complex_multiply(A_r[1], A_i[1], B_r[2], B_i[2]) - complex_multiply(A_r[2], A_i[2], B_r[1], B_i[1]),
)
end

"""
    Claculate cross product of two vectors. A and B represents complex vectors.
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product(A::AbstractMatrix{Real}, B::AbstractMatrix{Real})
    return complex_vector_cross_product(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent real and complex vectors, respectively.
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product(A::AbstractVector{Real}, B::AbstractMatrix{Real})
    return complex_vector_cross_product(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent complex and real vectors, respectively.
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product(A::AbstractMatrix{Real}, B::AbstractVector{Real})
    return complex_vector_cross_product(A[:,1], A[:,2], B, zero(B))
end

# Vector cross product, returns SMatrix ------------------------------------------------------------------------------------------
"""
    Claculate cross product of two vectors, inputs are real and imaginary parts of the two vectors. Returns SMatrix
Each is represented as 3-element array.
"""
function complex_vector_cross_product_SMatrix(A_r::AbstractVector{Real}, A_i::AbstractVector{Real}, B_r::AbstractVector{Real}, B_i::AbstractVector{Real})
    return vcat(
    complex_multiply_SMatrix(A_r[2], A_i[2], B_r[3], B_i[3]) - complex_multiply_SMatrix(A_r[3], A_i[3], B_r[2], B_i[2]),
    complex_multiply_SMatrix(A_r[3], A_i[3], B_r[1], B_i[1]) - complex_multiply_SMatrix(A_r[1], A_i[1], B_r[3], B_i[3]),
    complex_multiply_SMatrix(A_r[1], A_i[1], B_r[2], B_i[2]) - complex_multiply_SMatrix(A_r[2], A_i[2], B_r[1], B_i[1]),
)
end

"""
    Claculate cross product of two vectors. A and B represents complex vectors. Returns SMatrix
A represents a complex vector and is represented as 3x2 Matrix, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product_SMatrix(A::AbstractMatrix{Real}, B::AbstractMatrix{Real})
    return complex_vector_cross_product_SMatrix(A[:,1], A[:,2], B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent real and complex vectors, respectively. Returns SMatrix
A represents a real vector and is represented as 3-element vector, while B represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product_SMatrix(A::AbstractVector{Real}, B::AbstractMatrix{Real})
    return complex_vector_cross_product_SMatrix(A, zero(A), B[:,1], B[:,2])
end

"""
    Claculate cross product of two vectors. A and B represent complex and real vectors, respectively. Returns SMatrix
B represents a real vector and is represented as 3-element vector, while A represents a complex vector and is represented as 3x2 Matrix
"""
function complex_vector_cross_product_SMatrix(A::AbstractMatrix{Real}, B::AbstractVector{Real})
    return complex_vector_cross_product_SMatrix(A[:,1], A[:,2], B, zero(B))
end

# Complex matrix inverse ------------------------------------------------------------------------------------------
"""
    Calculate the inverse of a complex matrix A+iB, where A and B are real matrices

There are two possible versions:
Version #1: This one preallocate. Faster and use more memory
for 1M evaluations: 1.771474 seconds (14.18 M allocations: 4.510 GiB, 4.52% gc time, 6.40% compilation time)
`
function complex_matrix_inversion(A, B)
    inv_A_B = inv(A) * B
    C = inv(A + B * inv_A_B)
    D = -inv_A_B * C
    return C, D
end
`

Version #2: No preallocation. Slower and use less memory
for 1M evaluations: 2.330454 seconds (19.00 M allocations: 6.512 GiB, 2.71% gc time)
`
function complex_matrix_inversion(A, B)
    C = inv(A + B * inv(A) * B)
    D = -1 .* inv(A) * B * C
    return C, D
end
`
"""
function complex_matrix_inversion(A::AbstractMatrix{Real}, B::AbstractMatrix{Real})
    @fastmath @inbounds inv_A_B = inv(A) * B
    @fastmath @inbounds C = inv(A + B * inv_A_B)
    @fastmath @inbounds D = -inv_A_B * C
    return C, D
end

# Complex matrix multiplication ------------------------------------------------------------------------------------------
"""
    Multiply two matrices A+iB and C+iD, where A,B,C,D are real matrices
"""
function complex_matrix_multiplication(A::AbstractMatrix{Real}, B::AbstractMatrix{Real}, C::AbstractMatrix{Real}, D::AbstractMatrix{Real})
    return @fastmath @inbounds (A * C - B * D, A * D + B * C)
end

end # module


text:
using ITensors

function site_inds(ny, nx, d=1)
  sh = Matrix{Index}(undef, ny, nx)
  sv = Matrix{Index}(undef, ny, nx)
  for iy in 1:ny, ix in 1:nx
    sh[iy, ix] = Index(d, "site,horiz,x=$ix,y=$iy")
    sv[iy, ix] = Index(d, "site,vert,x=$ix,y=$iy")
  end
  return sh, sv
end

function link_inds(ny, nx, d=1)
  ll = Matrix{Index}(undef, ny, nx)
  lr = Matrix{Index}(undef, ny, nx)
  lu = Matrix{Index}(undef, ny, nx)
  ld = Matrix{Index}(undef, ny, nx)
  for iy in 1:ny, ix in 1:nx
    ll[iy, ix] = Index(d, "link,left,x=$ix,y=$iy")
    lr[iy, ix] = Index(d, "link,right,x=$ix,y=$iy")
    lu[iy, ix] = Index(d, "link,up,x=$ix,y=$iy")
    ld[iy, ix] = Index(d, "link,down,x=$ix,y=$iy")
  end
  return ll, lr, lu, ld
end

per(n, N) = mod(n - 1, N) + 1

function ctmrg_environment((sh, sv))
  ny, nx = size(sh)
  Clu = Matrix{ITensor}(undef, ny, nx)
  Cru = Matrix{ITensor}(undef, ny, nx)
  Cld = Matrix{ITensor}(undef, ny, nx)
  Crd = Matrix{ITensor}(undef, ny, nx)

  Al = Matrix{ITensor}(undef, ny, nx)
  Ar = Matrix{ITensor}(undef, ny, nx)
  Au = Matrix{ITensor}(undef, ny, nx)
  Ad = Matrix{ITensor}(undef, ny, nx)

  ll, lr, lu, ld = link_inds(ny, nx)
  for iy in 1:ny, ix in 1:nx
    Clu[iy, ix] = randomITensor(ll[iy, ix], lu[iy, ix])
    Cru[iy, ix] = randomITensor(lr[iy, ix], lu[iy, ix])
    Cld[iy, ix] = randomITensor(ll[iy, ix], ld[iy, ix])
    Crd[iy, ix] = randomITensor(lr[iy, ix], ld[iy, ix])
    iyp, ixp = per(iy + 1, ny), per(ix + 1, nx)
    Al[iy, ix] = randomITensor(sh[iy, ix], ll[iy, ix], ll[iyp, ix])
    Ar[iy, ix] = randomITensor(sh[iy, ix], lr[iy, ix], lr[iyp, ix])
    Au[iy, ix] = randomITensor(sv[iy, ix], lu[iy, ix], lu[iy, ixp])
    Ad[iy, ix] = randomITensor(sv[iy, ix], ld[iy, ix], ld[iy, ixp])
  end
  normalize!((Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad))
  return (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad)
end

function calc_κ(iy, ix, T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir="left")
  ny, nx = size(T)
  iyp, ixp = per(iy + 1, ny), per(ix + 1, nx)
  normC = scalar(Clu[iy, ix] * Cld[iy, ix] * Cru[iy, ix] * Crd[iy, ix])
  normAlr = scalar(
    Clu[iy, ix] * Cru[iy, ix] * Al[iy, ix] * Ar[iy, ix] * Cld[iyp, ix] * Crd[iyp, ix]
  )
  normAud = scalar(
    Clu[iy, ix] * Cld[iy, ix] * Au[iy, ix] * Ad[iy, ix] * Cru[iy, ixp] * Crd[iy, ixp]
  )
  normT = scalar(
    Clu[iy, ix] *
    Al[iy, ix] *
    Cld[iyp, ix] *
    Au[iy, ix] *
    T[iy, ix] *
    Ad[iyp, ix] *
    Cru[iy, ixp] *
    Ar[iy, ixp] *
    Crd[iyp, ixp],
  )
  return normT, normAlr, normAud, normC
end

function calc_κ(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir="left")
  ny, nx = size(T)
  κ = Matrix{Float64}(undef, ny, nx)
  for iy in 1:ny, ix in 1:nx
    normT, normAlr, normAud, normC = calc_κ(
      iy, ix, T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir=dir
    )
    κ[iy, ix] = normT * normC / (normAlr * normAud)
  end
  return κ
end

function normalize!((Clu, Cru, Cld, Crd))
  ny, nx = size(Clu)
  for iy in 1:ny, ix in 1:nx
    Clu[iy, ix] /= norm(Clu[iy, ix])
    Cld[iy, ix] /= norm(Cld[iy, ix])
    Cru[iy, ix] /= norm(Cru[iy, ix])
    Crd[iy, ix] /= norm(Crd[iy, ix])
    normC4 = scalar(Clu[iy, ix] * Cld[iy, ix] * Cru[iy, ix] * Crd[iy, ix])
    normC4 < 0 ? normClu = -abs(normC4)^(1 / 4) : normClu = normC4^(1 / 4)
    Clu[iy, ix] /= normClu
    Cld[iy, ix] /= abs(normClu)
    Cru[iy, ix] /= abs(normClu)
    Crd[iy, ix] /= abs(normClu)
  end
end

function normalize!((Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir="left")
  normalize!((Clu, Cru, Cld, Crd))
  ny, nx = size(Clu)
  for iy in 1:ny, ix in 1:nx
    Al[iy, ix] /= norm(Al[iy, ix])
    Ar[iy, ix] /= norm(Ar[iy, ix])
    Au[iy, ix] /= norm(Au[iy, ix])
    Ad[iy, ix] /= norm(Ad[iy, ix])
    iyp, ixp = per(iy + 1, ny), per(ix + 1, nx)
    normAlr = scalar(
      Clu[iy, ix] * Cru[iy, ix] * Al[iy, ix] * Ar[iy, ix] * Cld[iyp, ix] * Crd[iyp, ix]
    )
    normAlr < 0 ? normAl = -abs(normAlr)^(1 / 2) : normAl = normAlr^(1 / 2)
    Al[iy, ix] /= normAl
    Ar[iy, ix] /= abs(normAl)
    normAud = scalar(
      Clu[iy, ix] * Cld[iy, ix] * Au[iy, ix] * Ad[iy, ix] * Cru[iy, ixp] * Crd[iy, ixp]
    )
    normAud < 0 ? normAu = -abs(normAud)^(1 / 2) : normAu = normAud^(1 / 2)
    Au[iy, ix] /= normAu
    Ad[iy, ix] /= abs(normAu)
  end
end

function leftright_move!(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir="left", maxdim=5)
  ny, nx = size(T)
  P = Vector{ITensor}(undef, ny)
  P⁻ = Vector{ITensor}(undef, ny)
  if dir == "left" || dir == "up"
    xrange = 1:nx
  elseif dir == "right" || dir == "down"
    xrange = per.((nx - 1):-1:0, nx)
  end
  for ix in xrange
    ixm = per(ix - 1, nx)
    ixp = per(ix + 1, nx)
    ixpp = per(ix + 2, nx)
    for iy in 1:ny
      iym = per(iy - 1, ny)
      iyp = per(iy + 1, ny)

      Cu =
        Al[iym, ix] *
        Clu[iym, ix] *
        Au[iym, ix] *
        T[iym, ix] *
        Au[iym, ixp] *
        T[iym, ixp] *
        Cru[iym, ixpp] *
        Ar[iym, ixpp]
      @assert order(Cu) == 4
      Cd =
        Al[iy, ix] *
        Cld[iyp, ix] *
        Ad[iyp, ix] *
        T[iy, ix] *
        Ad[iyp, ixp] *
        T[iy, ixp] *
        Crd[iyp, ixpp] *
        Ar[iy, ixpp]
      @assert order(Cd) == 4
      if dir == "left" || dir == "up"
        li = commonindex(Cru[iy, ixpp], Crd[iy, ixpp])
        si = commonindex(Au[iy, ixp], Ad[iy, ixp])
      elseif dir == "right" || dir == "down"
        li = commonindex(Clu[iy, ix], Cld[iy, ix])
        si = commonindex(Au[iy, ix], Ad[iy, ix])
      end
      Cup = prime(Cu, (li, si))
      ρ = Cd * Cup
      if dir == "left" || dir == "right"
        utags = "$dir,link,x=$ixp,y=$iy"
      elseif dir == "up" || dir == "down"
        utags = "$dir,link,x=$iy,y=$ixp"
      end
      U, S, Vh, spec, u, v = svd(
        ρ, (li, si); utags=utags, vtags="tmp", maxdim=maxdim, cutoff=0.0
      )
      V = dag(Vh)
      U *= δ(u, v)
      invsqrtS = S
      for i in 1:dim(u)
        invsqrtS[i, i] = inv(sqrt(S[i, i]))
      end
      P[iy] = Cup * V * invsqrtS
      P⁻[iy] = Cd * dag(U) * invsqrtS
    end
    for iy in 1:ny
      iym = per(iy - 1, ny)
      iyp = per(iy + 1, ny)
      if dir == "left" || dir == "up"
        Al[iy, ixp] = Al[iy, ix] * P[iy] * T[iy, ix] * P⁻[iyp]
        Clu[iy, ixp] = Clu[iy, ix] * Au[iy, ix] * P⁻[iy]
        Cld[iy, ixp] = Cld[iy, ix] * Ad[iy, ix] * P[iy]
      elseif dir == "right" || dir == "down"
        Ar[iy, ixp] = Ar[iy, ixpp] * P[iy] * T[iy, ixp] * P⁻[iyp]
        Cru[iy, ixp] = Cru[iy, ixpp] * Au[iy, ixp] * P⁻[iy]
        Crd[iy, ixp] = Crd[iy, ixpp] * Ad[iy, ixp] * P[iy]
      end
    end
  end
  return normalize!((Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad))
end

function swapdiag(M)
  Mp = permutedims(M, [2, 1])
  ny, nx = size(Mp)
  for iy in 1:ny, ix in 1:nx
    Mp[iy, ix] = M[ix, iy]
  end
  return Mp
end
function rotate_environment(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad))
  return swapdiag(T), swapdiag.((Clu, Cld, Cru, Crd)), swapdiag.((Au, Ad, Al, Ar))
end

printdiv() = println("\n****************************************")
printstepdiv() = println("\n##################################################")

function sweepsdims(stepsizes::Vector{Int}, dims::Vector{Int})
  nstep = length(stepsizes)
  maxdims = zeros(Int, stepsizes[end])
  for i in 1:stepsizes[1]
    maxdims[i] = dims[1]
  end
  for j in 2:nstep
    for i in (stepsizes[j - 1] + 1):stepsizes[j]
      maxdims[i] = dims[j]
    end
  end
  return maxdims
end

function check_environment(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad))
  ny, nx = size(T)
  for iy in 1:ny
    for ix in 1:nx
      @assert order(Clu[iy, ix]) == 2
      @assert order(Cru[iy, ix]) == 2
      @assert order(Cld[iy, ix]) == 2
      @assert order(Crd[iy, ix]) == 2
      @assert order(Al[iy, ix]) == 3
      @assert order(Ar[iy, ix]) == 3
      @assert order(Au[iy, ix]) == 3
      @assert order(Ad[iy, ix]) == 3
      @assert order(T[iy, ix]) == 4
      @assert order(T[iy, ix]) == 4
      @assert order(T[iy, ix]) == 4
      @assert order(T[iy, ix]) == 4
      @assert length(commoninds(Clu[iy, ix], Cru[iy, ix])) == 1
      @assert length(commoninds(Clu[iy, ix], Cld[iy, ix])) == 1
      @assert length(commoninds(Cld[iy, ix], Crd[iy, ix])) == 1
      @assert length(commoninds(Cru[iy, ix], Crd[iy, ix])) == 1
    end
  end
end

function ctmrg(T::Matrix{ITensor}, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); verbose=false)
  ny, nx = size(T)

  verbose && println("Original:")
  verbose && @show calc_κ(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad))

  nstep = 1000
  maxdim = 10
  dirs = ["left", "up", "right", "down"]
  for ctmrg_step in 1:nstep
    verbose && printstepdiv()
    dir = dirs[per(ctmrg_step, length(dirs))]
    verbose && @show ctmrg_step, dir

    if dir == "left" || dir == "right"
      leftright_move!(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir=dir, maxdim=maxdim)
    elseif dir == "up" || dir == "down"
      T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad) = rotate_environment(
        T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad)
      )
      leftright_move!(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir=dir, maxdim=maxdim)
      T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad) = rotate_environment(
        T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad)
      )
    end

    check_environment(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad))

    verbose && @show Mκ = calc_κ(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir=dir)
    verbose && @show abs(prod(vec(Mκ)))^(1 / (nx * ny))
  end
  Mκ = calc_κ(T, (Clu, Cru, Cld, Crd), (Al, Ar, Au, Ad); dir=dir)
  κave = abs(prod(vec(Mκ)))^(1 / (nx * ny))
  return κave
end


text:
abstract RenderObject
typealias RenderStyle Dict{Symbol,Union(Integer,FloatingPoint,String)}

function kw_init(self::RenderObject, args...)
    for (key, value) in args2dict(args...)
        self.style[key] = value
    end
end

type LineObject <: RenderObject
    style::RenderStyle
    p::Point
    q::Point

    function LineObject(p, q, args...)
        self = new(RenderStyle(), p, q)
        kw_init(self, args...)
        self
    end
end

function boundingbox(self::LineObject, context)
    BoundingBox(self.p, self.q)
end

function draw(self::LineObject, context)
    line(context.draw, self.p.x, self.p.y, self.q.x, self.q.y)
end

type LabelsObject <: RenderObject
    style::RenderStyle
    points::Vector{Point}
    labels::Vector
    angle::Float64
    halign::ASCIIString
    valign::ASCIIString

    function LabelsObject(points, labels, args...; angle=0., halign="center", valign="center")
        self = new(RenderStyle(), points, labels, angle, halign, valign)
        kw_init(self, args...)
        self
    end
end

__halign_offset = [ "right"=>Vec2(-1,0), "center"=>Vec2(-.5,.5), "left"=>Vec2(0,1) ]
__valign_offset = [ "top"=>Vec2(-1,0), "center"=>Vec2(-.5,.5), "bottom"=>Vec2(0,1) ]

function boundingbox(self::LabelsObject, context)
    bb = BoundingBox()
    push_style(context, self.style)

    angle = self.angle * pi/180.

    height = textheight(context.draw, self.labels[1])
    ho = __halign_offset[self.halign]
    vo = __valign_offset[self.valign]

    for i = 1:length(self.labels)
        pos = self.points[i]
        width = textwidth(context.draw, self.labels[i])

        p = Point(pos.x + width * ho.x, pos.y + height * vo.x)
        q = Point(pos.x + width * ho.y, pos.y + height * vo.y)

        bb_label = BoundingBox(p, q)
        if angle != 0
            bb_label = rotate(bb_label, angle, pos)
        end
        bb += bb_label
    end

    pop_style(context)
    return bb
end

function draw(self::LabelsObject, context)
    for i in 1:length(self.labels)
        p = self.points[i]
        text(context.draw, p.x, p.y, self.labels[i];
             angle=self.angle, halign=self.halign, valign=self.valign)
    end
end

type CombObject <: RenderObject
    style::RenderStyle
    points::Vector{Point}
    dp

    function CombObject(points, dp, args...)
        self = new(RenderStyle())
        kw_init(self, args...)
        self.points = points
        self.dp = dp
        self
    end
end

function boundingbox(self::CombObject, context::PlotContext)
    return BoundingBox(self.points...)
end

function draw(self::CombObject, context::PlotContext)
    for p in self.points
        move_to(context.draw, p.x, p.y)
        rel_line_to(context.draw, self.dp.x, self.dp.y)
    end
    stroke(context.draw)
end

type SymbolObject <: RenderObject
    style::RenderStyle
    pos::Point

    function SymbolObject(pos, args...)
        self = new(RenderStyle(), pos)
        kw_init(self, args...)
        self
    end
end

function boundingbox(self::SymbolObject, context)
    push_style(context, self.style)
    symbolsize = get(context.draw, "symbolsize")
    pop_style(context)

    x = self.pos.x
    y = self.pos.y
    d = 0.5*symbolsize
    return BoundingBox(x-d, x+d, y-d, y+d)
end

function draw(self::SymbolObject, context)
    symbols(context.draw, [self.pos.x], [self.pos.y])
end

type SymbolsObject <: RenderObject
    style::RenderStyle
    x
    y

    function SymbolsObject(x, y, args...)
        self = new(RenderStyle())
        kw_init(self, args...)
        self.x = x
        self.y = y
        self
    end
end

function boundingbox(self::SymbolsObject, context::PlotContext)
    xmin = min(self.x)
    xmax = max(self.x)
    ymin = min(self.y)
    ymax = max(self.y)
    return BoundingBox(xmin, xmax, ymin, ymax)
end

function draw(self::SymbolsObject, context::PlotContext)
    symbols(context.draw, self.x, self.y)
end

type TextObject <: RenderObject
    style::RenderStyle
    pos::Point
    str::String
    angle::Float64
    halign::ASCIIString
    valign::ASCIIString

    function TextObject(pos, str, args...; angle=0., halign="center", valign="center")
        self = new(RenderStyle(), pos, str, angle, halign, valign)
        kw_init(self, args...)
        self
    end
end

function boundingbox(self::TextObject, context::PlotContext)
    push_style(context, self.style)
    angle = self.angle * pi/180.
    width = textwidth(context.draw, self.str)
    height = textheight(context.draw, self.str)
    pop_style(context)

    hvec = width * __halign_offset[self.halign]
    vvec = height * __valign_offset[self.valign]

    bb = BoundingBox(self.pos.x + hvec.x, self.pos.x + hvec.y,
                     self.pos.y + vvec.x, self.pos.y + vvec.y)
    bb = rotate(bb, angle, self.pos)
    return bb
end

function draw(self::TextObject, context::PlotContext)
    text(context.draw, self.pos.x, self.pos.y, self.str;
         angle=self.angle, halign=self.halign, valign=self.valign)
end

type PathObject <: RenderObject
    style::RenderStyle
    x::Vector{Float64}
    y::Vector{Float64}

    function PathObject(x, y, args...)
        self = new(RenderStyle())
        kw_init(self, args...)
        self.x = x
        self.y = y
        self
    end
end

function boundingbox(self::PathObject, context)
    xmin = min(self.x)
    xmax = max(self.x)
    ymin = min(self.y)
    ymax = max(self.y)
    return BoundingBox(xmin, xmax, ymin, ymax)
end

function draw(self::PathObject, context)
    curve(context.draw, self.x, self.y)
end

type PolygonObject <: RenderObject
    style::RenderStyle
    points::Vector{Point}

    function PolygonObject(points, args...)
        self = new(RenderStyle())
        kw_init(self, args...)
        self.points = points
        self
    end
end

function boundingbox(self::PolygonObject, context)
    return BoundingBox(self.points...)
end

function draw(self::PolygonObject, context)
    polygon(context.draw, self.points)
end

type ImageObject <: RenderObject
    style::RenderStyle
    img
    bbox

    function ImageObject(img, bbox, args...)
        self = new(RenderStyle(), img, bbox)
        kw_init(self, args...)
        self
    end
end

function boundingbox(self::ImageObject, context)
    return self.bbox
end

function draw(self::ImageObject, context)
    ll = lowerleft(self.bbox)
    w = width(self.bbox)
    h = height(self.bbox)
    if context.geom.yflipped
        image(context.draw, self.img, ll.x, ll.y+h, w, -h)
    else
        image(context.draw, self.img, ll.x, ll.y, w, h)
    end
end

# defaults

#function boundingbox(self::RenderObject, context)
#    return BoundingBox()
#end

function render(self::RenderObject, context)
    push_style(context, self.style)
    draw(self, context)
    pop_style(context)
end


text:
# This file has been generated automatically by `generate_compile_time_mpi_constants`.
# Do not edit.

# MPI library version (for reference only; this is not a constant)
const MPI_LIBRARY_VERSION_ = "Microsoft MPI i686"

# Compile-time constants

# Version:
# (We add an underscore to avoid a name conflict with `MPI.jl`)
const MPI_VERSION_ = Cint(2)
const MPI_SUBVERSION_ = Cint(0)

# Implementation limits:
const MPI_MAX_DATAREP_STRING = Cint(128)
const MPI_MAX_ERROR_STRING = Cint(512)
const MPI_MAX_INFO_KEY = Cint(255)
const MPI_MAX_INFO_VAL = Cint(1024)
const MPI_MAX_LIBRARY_VERSION_STRING = Cint(8192)
const MPI_MAX_OBJECT_NAME = Cint(128)
const MPI_MAX_PORT_NAME = Cint(256)
const MPI_MAX_PROCESSOR_NAME = Cint(128)

# Types

# Various (signed) integer types:
const MPI_Aint = Int32
const MPI_Count = Int64
const MPI_Fint = Int32
const MPI_Offset = Int64

# Status:
struct MPI_Status
    private0::Cint
    private1::Cint
    MPI_SOURCE::Cint
    MPI_TAG::Cint
    MPI_ERROR::Cint
end

# Opaque handles:
const MPI_Comm = UInt32
const MPI_Datatype = UInt32
const MPI_Errhandler = UInt32
const MPI_File = UInt32
const MPI_Group = UInt32
const MPI_Info = UInt32
const MPI_Message = UInt32
const MPI_Op = UInt32
const MPI_Request = UInt32
const MPI_Win = UInt32

# Function pointers:
const MPI_Comm_copy_attr_function = Ptr{Cvoid}
const MPI_Comm_delete_attr_function = Ptr{Cvoid}
const MPI_Comm_errhandler_function = Ptr{Cvoid}
const MPI_Comm_errhandler_fn = MPI_Comm_errhandler_function
const MPI_Copy_function = MPI_Comm_copy_attr_function
const MPI_Datarep_conversion_function = Ptr{Cvoid}
const MPI_Datarep_extent_function = Ptr{Cvoid}
const MPI_Delete_function = Ptr{Cvoid}
const MPI_File_errhandler_function = Ptr{Cvoid}
const MPI_File_errhandler_fn = Ptr{Cvoid}
const MPI_Grequest_cancel_function = Ptr{Cvoid}
const MPI_Grequest_free_function = Ptr{Cvoid}
const MPI_Grequest_query_function = Ptr{Cvoid}
const MPI_Type_copy_attr_function = Ptr{Cvoid}
const MPI_Type_delete_attr_function = Ptr{Cvoid}
const MPI_User_function = Ptr{Cvoid}
const MPI_Win_copy_attr_function = Ptr{Cvoid}
const MPI_Win_delete_attr_function = Ptr{Cvoid}
const MPI_Win_errhandler_function = Ptr{Cvoid}
const MPI_Win_errhandler_fn = MPI_Win_errhandler_function


text:
using FromFile
@from "Core.jl" import RecordType

"Assumes that `options` holds the user options::Options"
macro recorder(ex)
    quote
        if $(esc(:options)).recorder
            $(esc(ex))
        end
    end
end

function find_iteration_from_record(key::String, record::RecordType)
    iteration = 0
    while haskey(record[key], "iteration$(iteration)")
        iteration += 1
    end
    return iteration - 1
end


text:
using LearningStrategy
using DataFrames

import LearningStrategy: hook

struct DFTracer <: LearningStrategy
  f
  b::Int
  df::DataFrame

  DFTracer(T::Vector{<:Type}, names::Vector{Symbol}, f, b = 1) =
    new(f, b, DataFrame(T, names, 0))
end

function hook(t::DFTracer, m, i)
  (mod1(i, t.b) == t.b) && push!(t.df, t.f(m, i))
  nothing
end

dftracer = DFTracer([Float64, Symbol], [:reward, :action], (x, i) -> begin
  # ...
  (42, :magic)
end)


text:
module GMDApp

import Base: parse, getindex
using DelimitedFiles: readdlm

using Gaston

using Statistics

const DEFAULT_INPUT = "input.txt"
const WARNING_INPUT_1 =  "
WARNING: $(DEFAULT_INPUT) found in the GMDApp home directory.
Default parameters will be used during the execution.
Modify or delete $(DEFAULT_INPUT) from the GMDApp home directory
for user-defined input.
"
const WARNING_INPUT_2 =  "
WARNING: $(DEFAULT_INPUT) not found in the GMDApp home directory.
Introduce the name of your input file: "
const WELCOME_MESSAGE = "
GMDApp -- Ground Motion Database Search Application
Version 0.2.0 64-Bit
Copyright (c) 2022 Diego Casas, Gustavo Araujo, and Alexander Arciniegas.
MIT License.
This is a short implementation in Julia of the PEER Ground Motion Database
Web Application Algorithm for Time Series Selection.
The app is an extension of the PEER procedure for user-defined databases.
"
struct Record
    RSN::Int64
    EQID::Int64
    EQName::String
    Year::Int64
    SSN::Int64
    SName::String
    FaultType::Int64
    Magnitude::Float64
    Rjb::Float64
    ClstD::Float64
    Vs30::Float64
    PGA::Float64
    PGV::Float64
    PGD::Float64
    Dur::Float64
    S::Vector{Float64}
    FileNameHorizontal1::String
    FileNameHorizontal2::String
    FileNameVertical::String
    function Record(;
        RSN,
        EQID,
        EQName,
        Year,
        SSN,
        SName,
        FaultType,
        Magnitude,
        Rjb,
        ClstD,
        Vs30,
        PGA,
        PGV,
        PGD,
        Dur,
        S,
        FileNameHorizontal1,
        FileNameHorizontal2,
        FileNameVertical,
    )
        new(
            RSN,
            EQID,
            EQName,
            Year,
            SSN,
            SName,
            FaultType,
            Magnitude,
            Rjb,
            ClstD,
            Vs30,
            PGA,
            PGV,
            PGD,
            Dur,
            S,
            FileNameHorizontal1,
            FileNameHorizontal2,
            FileNameVertical,
        )
    end
end

function getindex(rec::Record, k::AbstractString)
    getproperty(rec, Symbol(k))
end

function parse(::Type{Tuple{T}}, s::AbstractString) where {T}
    Tuple(parse.(T, split(s, ',')))
end

"""
    read_header(hline, nt)

Returns a dictionary with field name-column number pairs.
Currently, the output dictionary only works when the records
file has field names equal to the Record struct fields.
"""
function read_header(hline::AbstractString, nt::Int64)
    fields = Dict()
    header = split(hline, ',')
    n = length(header)
    push!(fields, "FileNameVertical" => n)
    push!(fields, "FileNameHorizontal2" => n - 1)
    push!(fields, "FileNameHorizontal1" => n - 2)
    push!(fields, "S" => (n-nt-2):(n-3))
    for i in 1:(n-nt-3)
        push!(fields, header[i] => i)
    end
    fields
end

"""
    read_records(filename)

Returns a vector of Record instances read from the given file.
Also returns the periods vector for the spectra in the records.
"""
function read_records(filename)
    open(filename) do io
        nrec = countlines(io) - 3
        seekstart(io)
        readline(io) # title
        T = parse.(Float64, filter((~) ∘ isempty, split(readline(io), ',')))
        fields = read_header(readline(io), length(T))
        records = Vector{Record}(undef, nrec)
		for (i, row) in enumerate(eachrow(readdlm(io, ',', String)))
        # for (i, line) in enumerate(eachline(io))
            # row = split(line, ',')
            records[i] = Record(
                RSN = parse(Int64, row[fields["RSN"]]),
                EQID = parse(Int64, row[fields["EQID"]]),
                EQName = row[fields["Earthquake Name"]],
                Year = parse(Int64, row[fields["Year"]]),
                SSN = parse(Int64, row[fields["SSN"]]),
                SName = row[fields["Station Name"]],
                FaultType = parse(Int64, row[fields["FaultType"]]),
                Magnitude = parse(Float64, row[fields["Magnitude"]]),
                Rjb = parse(Float64, row[fields["Rjb"]]),
                ClstD = parse(Float64, row[fields["ClstD"]]),
                Vs30 = parse(Float64, row[fields["Vs30"]]),
                PGA = parse(Float64, row[fields["PGA"]]),
                PGV = parse(Float64, row[fields["PGV"]]),
                PGD = parse(Float64, row[fields["PGD"]]),
                Dur = parse(Float64, row[fields["Dur"]]),
                S = parse.(Float64, row[fields["S"]]),
                FileNameHorizontal1 = row[fields["FileNameHorizontal1"]],
                FileNameHorizontal2 = row[fields["FileNameHorizontal2"]],
                FileNameVertical = row[fields["FileNameVertical"]],
            )
        end
        records, T
    end
end

"""
    read_filters(filename)

Returns a function that is suitable for filtering an array of Record
instances from an input filter file. The output function checks for the
Record to satisfy all the conditions given in the input filter file.
It also check for all the values of the spectrum to be positive.
"""
function read_filters(filename)
    filters = Dict()
    postfilters = Dict()
    open(filename) do io
        for line in eachline(io)
            k, v = split(line)
            if k == "Nmax"
                Nmax = parse(Int64, v)
                push!(postfilters, "Nmax" => Nmax)
                continue
            end
            if k == "ScaleFactor"
                ScaleFactor = v == "NA" ? nothing : parse(Tuple{Float64}, v)
                push!(postfilters, "ScaleFactor" => ScaleFactor)
                continue
            end
            if v != "NA"
                T = k == "FaultType" ? Int64 : Float64
                push!(filters, k => parse(Tuple{T}, v))
            end
        end
    end
    function f(rec)
        for (k, v) in filters
            # returns false at the firt false condition
            # i.e. ALL conditions must be satified
            if k == "FaultType"
                !(rec[k] in v) && return false
            else
                !(first(v) ≤ rec[k] ≤ last(v)) && return false
            end
        end
        # check that all the spectrum is positive
        !any(rec.S .≤ 0)
    end
    f, postfilters
end

"""
    read_target_spectrum(filename)

Returns the periods and spectral acceleration of the target spectrum (T, S),
and arrays for the application for weights (T_w, w).
"""
function read_target_spectrum(filename)
    open(filename) do io
        readline(io) # title
        T_w = parse.(Float64, filter((~) ∘ isempty, split(readline(io), ',')))
        w = parse.(Float64, filter((~) ∘ isempty, split(readline(io), ',')))
        data = readdlm(io, ',', String, skipstart = 1)
        T = parse.(Float64, data[:, 1])
        S = parse.(Float64, data[:, 2])
        T, S, T_w, w
    end
end

function loginterp(x, xp, yp)
    n = length(x)
    m = length(xp)
    y = Vector{Float64}(undef, n)
    j = 1
    first_xp = first(xp)
    first_yp = first(yp)
    last_xp = last(xp)
    last_yp = last(yp)
    for i in 1:n
        if x[i] < first_xp
            y[i] = first_yp
        elseif x[i] > last_xp
            y[i] = last_yp
        else
            while j < m
                if xp[j] ≤ x[i] ≤ xp[j+1]
                    k = log(yp[j+1] / yp[j]) / log(xp[j+1] / xp[j])
                    y[i] = exp(k * log(x[i] / xp[j]) + log(yp[j]))
                    # The next point could be in the same bracket,
                    # so we break the loop without accumulating j.
                    break
                end
                j += 1
            end
        end
    end
    y
end

function read_input()
    files = Dict()
    if isfile(DEFAULT_INPUT)
        filename = DEFAULT_INPUT
        print(WARNING_INPUT_1)
        println()
    else
        print(WARNING_INPUT_2)
        filename = readline()
        println()
    end
    open(filename) do io
        for line in eachline(io)
            k, v = split(line)
            push!(files, k => v)
        end
    end
    files
end

function genplot(T, S_target, T_target, S_best, S_mean, S_mean_plus_std, S_mean_minus_std)
    n = length(T)
    nrec = size(S_best)[2]
    plt = plot(T, S_best[:, 1], w = "l", legend =  "''", lw = 2, lc = :gray,
        Axes(grid = :on,
          key = "left",
          axis = "loglog",
          xrange = (T_target[1], T_target[end]),
          xlabel = "'T [s]'",
          ylabel = "'Sa [g]'",
          )
    )
    for i = 1:nrec
        plot!(T, S_best[:, i], w = "l", legend =  "''", lw = 2, lc = :gray)
    end
    plot!(T, S_target, w = "l", legend = "'Target'",  lw = 2, lc = :red)
    plot!(T, S_mean, w = "l", legend = "'Mean'",  lw = 2, lc = :black)
    plot!(T, S_mean_plus_std, w = "l", legend = "'Mean+Std'",  dt = 2, lw = 2, lc = :black)
    plot!(T, S_mean_minus_std,  w = "l", legend = "'Mean-Std'",  dt = 2, lw = 2, lc = :black)
    plt
end

function write_criteria(io, filterfile, T_w, w_orig)
    println(io, "-- Summary of GMDApp Search Criteria --")
    open(filterfile) do io2
        for line in eachline(io2)
            k, v = split(line)
            println(io, k, ":,", v)
        end
    end
    print(io, "Period Array:,")
    for v in T_w
        print(io, v, ",")
    end
    println(io)
    print(io, "Weight Array:,")
    for v in w_orig
        print(io, v, ",")
    end
    println(io)
end

function write_summary(io, records, e², f)
    println(io, "-- Summary of Metadata of Selected Records --")
    if io == stdout
        println(
            io,
            "ID, RSN, EQID, MSE, ScaleFactor, Dur, Magnitude, FaultType, Rjb, ClstD, Vs30",
        )
        for j in 1:length(records)
            rec = records[j]
            println(
                io, j, ", " ,rec.RSN, ", ", rec.EQID, ", ", round(e²[j], sigdigits = 4), ", ",
                round(f[j], sigdigits = 2), ", ", rec.Dur, ", ",  rec.Magnitude,
                ", ", rec.FaultType, ", ", rec.Rjb, ", ", rec.ClstD, ", ",
                rec.Vs30,
            )
        end
    else
        println(
            io,
            "ID,RSN,MSE,ScaleFactor,Dur,Earthquake Name,Year,Station Name,Magnitude,FaultType,Rjb,ClstD,Vs30,FileNameHorizontal1,FileNameHorizontal2,FileNameVertical",
        )
        for j in 1:length(records)
            rec = records[j]
            println(
                io, j, "," ,rec.RSN, ",", e²[j], ",", f[j], ",", rec.Dur, ",",
                rec.EQName, ",", rec.Year, ",", rec.SName, ",",
                rec.Magnitude, ",", rec.FaultType, ",", rec.Rjb, ",",
                rec.ClstD, ",", rec.Vs30, ",", rec.FileNameHorizontal1, ",",
                rec.FileNameHorizontal2, ",", rec.FileNameVertical,
            )
        end
    end
end

function write_scaled_spectra(
    io,
    records,
    T,
    S_target,
    S_mean,
    S_mean_plus_std,
    S_mean_minus_std,
    S_best,
)
    println(io, "-- Scaled Spectra used in Search & Scaling --")
    # write header
    print(io, "Period,TargetSpectrum,Mean,Mean+StdDev,Mean-StdDev")
    for rec in records
        print(io, ",RSN-$(rec.RSN)")
    end
    println(io)
    # write lines with data
    for i in 1:length(T)
        print(
            io,
            T[i],
            ",",
            S_target[i],
            ",",
            S_mean[i],
            ",",
            S_mean_plus_std[i],
            ",",
            S_mean_minus_std[i],
        )
        for j in 1:length(records)
            print(io, ",", S_best[i, j])
        end
        println(io)
    end
end

function write_unscaled_spectra(io, records, T)
    println(io, "-- Unscaled Spectra --")
    # write header
    print(io, "Period")
    for rec in records
        print(io, ",RSN-$(rec.RSN)")
    end
    println(io)
    # write lines with data
    for i in 1:length(T)
        print(io, T[i])
        for rec in records
            print(io, ",", rec.S[i])
        end
        println(io)
    end
end

function main()
    print(WELCOME_MESSAGE)
    println()
    files = read_input()

    records, T = read_records(files["flatfile"])
    n = length(T)

    criteria, postfilters = read_filters(files["filters"])
    filter!(criteria, records)
    nrec = length(records)
    if nrec == 0
        println("None of the records satisfy the given filters.")
        return
    end

    T_target, S_target, T_w, w_orig = read_target_spectrum(files["spectrum"])
    S_target = loginterp(T, T_target, S_target)
    w = loginterp(T, T_w, w_orig)
    w[(T.<first(T_target)).&(T.>last(T_target))] .= 0
    w ./= sum(w) # normalize weights

    f = Vector{Float64}(undef, nrec)
    e² = Vector{Float64}(undef, nrec)
    fmin = 0
    fmax = 9999
    ScaleFactor = postfilters["ScaleFactor"]
    if !isnothing(ScaleFactor)
        fmin = first(ScaleFactor)
        fmax = last(ScaleFactor)
    end
    for (i, rec) in enumerate(records)
        fi = exp(sum(w .* log.(S_target ./ rec.S)))
        if fi < fmin
            fi = fmin
        end
        fi < fmin && (fi = fmin)
        fi > fmax && (fi = fmax)
        f[i] = fi
        e²[i] = sum(w .* log.(S_target ./ (f[i] .* rec.S)) .^ 2)
    end

    k = postfilters["Nmax"]

    nrec = length(records)
    if nrec == 0
        println("None of the records satisfy the given Scale Factor range.")
        return
    end
    k = min(k, nrec)

    ibest = partialsortperm(e², 1:k)
    write_criteria(stdout, files["filters"], T_w, w_orig)
    write_summary(stdout, records[ibest], e²[ibest], f[ibest])

    # generate output
    S_best = Matrix{Float64}(undef, n, length(ibest))
    for (j, i) in enumerate(ibest)
        rec = records[i]
        S_best[:, j] .= f[i] .* rec.S
    end
    S_mean = Vector{Float64}(undef, n)
    S_std = Vector{Float64}(undef, n)
    for i in 1:n
        S_mean[i] = mean(S_best[i, :])
        S_std[i] = exp(std(log.(S_best[i, :])))
    end
    S_mean_plus_std = S_mean .* S_std
    S_mean_minus_std = S_mean ./ S_std
    plt = genplot(T, S_target, T_target, S_best, S_mean, S_mean_plus_std, S_mean_minus_std)

    open(files["output"], "w") do io
        write_criteria(io, files["filters"], T_w, w_orig)
        println(io)
        write_summary(io, records[ibest], e²[ibest], f[ibest])
        println(io)
        write_scaled_spectra(
            io,
            records[ibest],
            T,
            S_target,
            S_mean,
            S_mean_plus_std,
            S_mean_minus_std,
            S_best,
        )
        println(io)
        write_unscaled_spectra(io, records[ibest], T)
    end
    display(plt)
    println("Press Enter to exit...")
    readline(stdin)
end

function julia_main()::Cint
    try
        main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

end # module

max_stars_repo_name: gaaraujo/GMDApp

max_issues_repo_name: gaaraujo/GMDApp

max_forks_repo_name: gaaraujo/GMDApp


text:
const hs_wall_rebound = [
    v -> [-v[1],  v[2]],
    v -> [ v[1], -v[2]]
]

function step!(hs :: HardSpheres, pc::ParticleCollisions, bc::BoxCollisions)

    ϵ = hs.ϵ
    n = hs.n

    tempsp, i1, i2 = dt_min_position(pc)
    tempsm, j1, j2 = dt_min_position(bc)

    if tempsp < tempsm

        for i in 1:n
            hs.q[i] = hs.q[i] .+ tempsp .* hs.v[i]
        end

        δv = hs.v[i2] - hs.v[i1]
        δx = hs.q[i2] - hs.q[i1]

        J = (δv'δx) ./ 2ϵ

        hs.v[i1] = hs.v[i1] .+ J .* δx ./ 2ϵ
        hs.v[i2] = hs.v[i2] .- J .* δx ./ 2ϵ

        pc.dt .-= tempsp
        reset!(pc, i1)
        reset!(pc, i2)

        bc.dt .-= tempsp
        reset!(bc, i1)
        reset!(bc, i2)

        compute_dt!(pc, i1, hs)
        compute_dt!(pc, i2, hs)

        t, i = compute_dt(hs, i1)
        bc.dt[i1, i] = t

        t, i = compute_dt(hs, i2)
        bc.dt[i2, i] = t

    else

        for i in 1:n
            hs.q[i] = hs.q[i] + tempsm * hs.v[i]
        end

        hs.v[j1] = hs_wall_rebound[j2](hs.v[j1])

        pc.dt .-= tempsm
        reset!(pc, j1)
        compute_dt!(pc, j1, hs)

        bc.dt .-= tempsm
        reset!(bc, j1)
        t, i = compute_dt(hs, j1)
        bc.dt[j1, i] = t

    end

    return min(tempsp, tempsm)

end


text:
@testset "kernel_2" begin
  using LightGraphs
  kg = kronecker_generator(10,2)
  g = kernel_1(kg)
  keys_ = @inferred key_sampling(g)
  @test length(keys_) <= 64
  kg = kronecker_generator(2,2)
  g = kernel_1(kg)
  keys_ = @inferred key_sampling(g)
  @test length(keys_) <= 64
  T = eltype(g)
  parent_  = @inferred  kernel_2(g,0x01)
  @test length(parent_) <= 2^2
end


text:
#!/usr/bin/env julia

using Luxor

using Test

function test_lux_bb(;center = O)
    bx = BoundingBox(box(O + center, 200, 200, :none))
    box(bx, :stroke)
    for i in 0:14
        pt = if i < 4
            sethue("yellow")
            ang = i * π / 2
            O + 40 .* (cos(ang), sin(ang))
        else
            sethue("grey50")
            randompoint((1.5bx)...)
        end
        arrow(O, pt)
        pt2 = pointcrossesboundingbox(pt, bx)
        circle(pt2, 3, :stroke)
    end
    snapshot()
end

function test_pointcrossesboundingbox(fname)
    Drawing(800,400,:rec)
    background("brown")
    origin(200,200)
    test_lux_bb()
    translate(400,0)
    test_lux_bb(;center = (-45, -25))
    snapshot(;fname=fname)
    println("...finished point crosses bounding box test, saved in $(fname)")
end

fname = "point-crosses-bounding-box-test.png"

test_pointcrossesboundingbox(fname)


text:
include("quadrature.jl")
include("Derivatives.jl")
include("TDist.jl")
include("transforms.jl")

"""
    createTensorGrid(example, meshtheta, meshlambda, type)

Define a function ``f`` from ``R^k`` to ``Mat(n, n)``, such that ``f(z_0)_{ij} = p(z_0|z, θ_i, λ_j)``,
where ``i`` and ``j`` range over the meshgrids over ``θ`` and ``λ``. Optional arg ``type`` is ""Gaussian""
by default. If ``type`` is "Turan", then use Gauss-Turan quadrature to integrate out ``0`` variable.
"""
function getTensorGrid(train::trainingData{T1, T2}, test::testingData{T1}, priorθ, priorλ, nodesWeightsθ::nodesWeights, nodesWeightsλ::nodesWeights, transform, quadtype::String) where T1 <:Array{Float64} where T2<:Array{Float64}
    nodesθ = nodesWeightsθ.nodes
    nodesλ = nodesWeightsλ.nodes
    weightsθ = nodesWeightsθ.weights
    weightsλ = nodesWeightsλ.weights
    X = train.X; X0 = test.X0; z = train.z; n = size(X, 1); p = size(X, 2)
    l1 = length(nodesθ); l2 = length(nodesλ); l3 = size(nodesWeightsθ.weights, 2)
    function func_fixed(θ::Float64)
        return funcθ(θ, train, test, quadtype)
    end
    #precompute buffer for θ-values
    theta_param_list = Array{Union{θ_params{Array{Float64, 2}, Cholesky{Float64,Array{Float64, 2}}}, θ_param_derivs{Array{Float64, 2}, Cholesky{Float64,Array{Float64, 2}}}}}(undef, l1)
    for i=1:l1
        theta_param_list[i] = func_fixed(nodesθ[i])
    end

    #initialize tensor grid of weights
    n1 = length(nodesθ); n2 = length(nodesλ); n3 = size(weightsθ, 2);
    weightsTensorGrid = zeros(n1, n2, n3) #tensor grid of weights
    for i = 1:n1
        for j = 1:n2
            funcs = posterior_theta(nodesθ[i], nodesλ[j], priorθ, priorλ, train, test, boxCoxObj, theta_param_list[i], quadtype)
            for k = 1:n3
                weightsTensorGrid[i, j, k] = weightsθ[i, k]*weightsλ[j]*funcs[k]
            end
        end
    end
    #compute matrix of Jacobians stably (using shift and scale)
    jacvals = zeros(1, n2)
    df = transform.df
    for i = 1:n2
        jacvals[i] = sum(log.(abs.(map( x-> df(x, nodesλ[i]), z))))
    end
    maxz = maximum(jacvals)
    for i = 1:n2
        jacvals[i] = jacvals[i]-maxz
    end
    jacTensorGrid = repeat((exp.(jacvals)).^(1-p/n), n1, 1, n3)

    weightsTensorGrid = weightsTensorGrid .* jacTensorGrid
    weightsTensorGrid =  weightsTensorGrid/sum(weightsTensorGrid)

    #tensor grid for PDF p(z0|theta, lambda, z)
    tgridpdf = Array{Any, 3}(undef, l1, l2, l3)
    #tensor grid for CDF P(z0|theta, lambda, z)
    tgridcdf = Array{Any, 3}(undef, l1, l2, l3)

    for i = 1:l1
        for j = 1:l2
            funcs = partial_theta(nodesθ[i], nodesλ[j], train, test, boxCoxObj, theta_param_list[i], quadtype)
            for k = 1:l3
                #tgrid[i, j, k] = funcs[k]
                tgridpdf[i, j, k] = funcs[1][k]
                tgridcdf[i, j, k] = funcs[2][k]
            end
        end
    end
    function evalTgrid_pdf(z0)
        res = Array{Float64, 3}(undef, l1, l2, l3)
        for i=1:l1
            for j = 1:l2
                for k =1:l3
                    res[i, j, k] = tgridpdf[i, j, k](z0)
                end
            end
        end
        return res
    end
    function evalTgrid_cdf(z0)
        res = Array{Float64, 3}(undef, l1, l2, l3)
        for i=1:l1
            for j = 1:l2
                for k =1:l3
                    res[i, j, k] = tgridcdf[i, j, k](z0)
                end
            end
        end
        return res
    end
    #product rule loop for Turan quadrature (k>1). In the k=1 case, this reduces to taking the dot
    # product between the tensor grid and weights grid.
    function pdf(z0)
        grid = evalTgrid_pdf(z0)
        res = 0
        for i = 1:n3
            for j = 1:i
                res = binomial(i, j) * dot(grid[:, :, j],  weightsTensorGrid[:, :, n3-j+1])
            end
        end
        return res
    end
    function cdf(z0)
        grid = evalTgrid_cdf(z0)
        res = 0
        for i = 1:n3
            for j = 1:i
                res = binomial(i, j) * dot(grid[:, :, j], weightsTensorGrid[:, :, n3-j+1])
            end
        end
        return res
    end
    return (pdf, cdf)
    #return z0 -> dot(evalTgrid_pdf(z0), weightsTensorGrid), z0 -> dot(evalTgrid_cdf(z0), weightsTensorGrid)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@everywhere function calcnan(xa::Float64, p::Int64, start = 0.0, fin = 4.0)
  λarr = collect(linspace(start, fin, p))
  const pconst = 2720 # calculated with findeq()
  const iterconst = 50000
  const ϵ = 1e-6
  yarr = SharedArray(Float64, (p, pconst))
  fill!(yarr, NaN)
  @sync @parallel for i = 1:p
    λ = λarr[i]
    x = xa
    for j = 1:iterconst
      x = λ*x*(1-x)
    end
    yarr[i, 1] = x # add first value
    m = 1
    found = false
    while !found # while we still have new values
      x = λ*x*(1-x) # iterate once more
      for k = 1:m # iterate through all found values
        if almostequal(x, yarr[i, k], ϵ) # and check if we already have the new one
          found = true
          break
        end
      end
      if !found
        m += 1
        yarr[i, m] = x # if not, add it to the array
      end
    end
  end
  return λarr, yarr
end


text:
using Plots, Serialization, LaTeXStrings, Printf

push!(LOAD_PATH,joinpath(pwd(),"src"));

using GP
using ModelicaPipe_GPSSM

###########################
# Main loop
###########################
# parameters [50-100-1000-100 funkar bra]
nx = 50 # number of segments of simulated pipe
n_pre = 100 # number of points for initial dataset
n_large = 1000 # number of points for large dataset
m_sparse = 50 # number of points for sparse dataset
Θ₀ = GP.hyperparameters([1.0, 1.0], [1.0, 1.0, 1.0], 0.01, [3,1,2])

# generate the GP
GP_pipe, ixl, ixm = ModelicaPipe_GPSSM.generate_GP(n_pre, n_large, m_sparse, Θ₀)
k_pipe = ModelicaPipe_GPSSM.k_pipe

# get validation data
Xv, yv, ixv, data_df, input_df, validation_df = ModelicaPipe_GPSSM.read_csv("val", 4:6, 2)
T_val = Array(validation_df[:,2:51])' .- 273.15

# create ranges for simulation
t_range = 1:1000 #size(input_df,1)-1

# number of simulations for monte carlo type simulation
nmc=1

# simulation with sparse GP
T_sim_mc, V_sim_mc = @time ModelicaPipe_GPSSM.forward_simulation(t_range, input_df, GP_pipe, k_pipe, nx, true, nmc)

# plots
gr(); save_results=false

# results at pipe outlet
Δt = 5; Δx = 10; x=50;
ModelicaPipe_GPSSM.plot_results_t(T_sim_mc, V_sim_mc, T_val, Δt, x, t_range, "Time [s]", "Temperature [C]", "plots/posterior_plot_t.pdf", save_results, false);

# monte carlo simulation results
ModelicaPipe_GPSSM.plot_results_t(T_sim_mc, V_sim_mc, T_val, Δt, x, t_range, "Time [s]", "Temperature [C]", "plots/posterior_plot_t_mc.pdf", save_results, true);

# spatial distribution
x_range=1:50; t=20;
ModelicaPipe_GPSSM.plot_results_x(T_sim_mc, V_sim_mc, T_val, Δx, x_range, t, "Position [m]", "Temperature [C]", "plots/posterior_plot_x.pdf", save_results);

max_stars_repo_name: smnssn/GPDH

max_issues_repo_name: smnssn/GPDH

max_forks_repo_name: smnssn/GPDH


text:
function load(f::File{format"OMETIFF"}; dropunused=true, inmemory=true)
    open(f) do s
        ret = load(s; dropunused=dropunused, inmemory=inmemory)
    end
end

"""
    load(io; dropunused, verbose, inmemory) -> ImageMetadata.ImageMeta

Load an OMETIFF file using the stream `io`.

**Arguments**
- `dropunused::Bool`: controls whether dimensions of length 1 are dropped
  automatically (default) or not.
- `verbose::Bool`: if true then prints a progress bar during loading
- `inmemory::Bool`: controls whether arrays are fully loaded into memory
  (default) or left on disk and specific parts only loaded when accessed.

!!! tip
    The `inmemory=false` flag currently returns a read-only view of the data on
    the disk for data integrity reasons. In order to modify the contents, you
    must copy the data into an in-memory container--at least until
    [#52](https://github.com/tlnagy/OMETIFF.jl/issues/52) is fixed--like so:

    ```
    copy(arr)
    ```
"""
function load(io::Stream{format"OMETIFF"}; dropunused=true, verbose = true, inmemory=true)
    if io.filename !== nothing && !occursin(".ome.tif", io.filename)
        throw(FileIO.LoaderError("OMETIFF", "Not an OME TIFF file!", ErrorException("")))
    end

    orig_file = read(io, TiffFile)
    summary = load_comments(orig_file)

    # load master OME-XML that contains all information about this dataset
    omexml = load_master_xml(orig_file)

    # find all images in this dataset, can either have the Image or Pixel tag
    containers = findall("//*[@DimensionOrder]", omexml, ["ns" => namespace(omexml)])

    pos_names = nodecontent.(findall("/ns:OME/ns:Image/ns:StageLabel[@Name]/@Name", omexml, ["ns"=>namespace(omexml)]))
    # if all position names aren't unique then substitute names
    if length(pos_names) == 0 || !allunique(pos_names)
        pos_names = ["Pos$i" for i in 1:length(containers)]
    end

    dimlist = []
    axeslist = []

    for (idx, container) in enumerate(containers)
        dims, axes_info = build_axes(container)
        push!(dimlist, dims)
        push!(axeslist, axes_info)
    end

    ifd_indices = OrderedDict{Int, NTuple{4, Int}}()
    ifd_files = OrderedDict{Int, Tuple{UUID, String}}()
    obs_filepaths = Dict{String, Int}()
    for (idx, container) in enumerate(containers)
        ifdindex!(ifd_indices, ifd_files, obs_filepaths, container, dimlist[idx], orig_file, idx)
    end

    files, ifds = get_ifds(orig_file, ifd_indices, ifd_files)

    # determine size of all the dims from the ifds in the tiff file instead
    # of the sizes embedded in the Pixel node
    true_dims = [Set{Int}() for i in 1:4]
    for ifd_index in keys(ifds), dim in 1:4
        push!(true_dims[dim], ifd_index[dim])
    end

    # generate new master dim list with the true dims from above
    master_dims = dimlist[1]
    dimnames = keys(master_dims)[3:6]
    new_data = [dimnames[i]=>length(true_dims[i]) for i in 1:4]
    master_dims = merge(master_dims, new_data)

    masteraxis = copy(axeslist[1])
    masteraxis[6] = Axis{:position}(to_symbol.(pos_names))

    # check if the axes computed earlier are the correct length, if not,
    # we're not sure about the information through the whole movie so lets
    # strip the units
    err_axes = findall(length.(masteraxis) .!== values(master_dims))
    for ax in err_axes
        masteraxis[ax] = masteraxis[ax](1:master_dims[ax])
    end

    elapsed_times = get_elapsed_times(containers, master_dims, masteraxis)

    if inmemory
        loaded = inmemoryarray(ifds, master_dims; verbose = verbose)
    else
        loaded = DiskOMETaggedImage(ifds, values(master_dims));
    end

    data = fixcolors(loaded, first(values(ifds))[2])

    # find dimensions of length 1 and remove them
    if dropunused
        unused_dims = findall(values(master_dims) .== 1)
        data = dropdims(data; dims=tuple(unused_dims...))
        deleteat!(masteraxis, unused_dims)
    end

    ImageMeta(AxisArray(data, masteraxis...),
              Description=summary,
              Elapsed_Times=elapsed_times)
end

"""
    inmemoryarray(ifds, dims, rawtype, mappedtype) -> Array

Builds an in-memory high-dimensional image using the mapping provided by `ifds`
from indices to [`OMETIFF.IFD`](@ref) objects. The IFD objects store handles to
the file objects and the offsets for the data. `dims` stores the size of each
named dimension. The `rawtype` parameter describes the storage layout of each
element on disk and `mappedtype` is the corresponding fixed or floating point
type.
"""
function inmemoryarray(ifds, dims::NamedTuple; verbose = true)

    ifd = first(values(ifds))[2]
    cache = getcache(ifd)

    data = similar(cache, dims...)

    freq = verbose ? 1 : Inf

    # iterate over each IFD
    @showprogress for (indices, (file, ifd)) in ifds
        read!(cache, file, ifd)
        data[:, :, indices...] .= cache'
    end

    return data
end