text:
module DataScripts

using LsqFit
using Plots
using DataFrames
using SFGTools

export plotrecent, pump_probe_process, cal_spectrometer

const DATA_PATH = "/Users/lackner/Documents/DataSFG"

# f = figure()
# close(f)

function plotrecent(idx=1; pixel=false, sfwl=false, dograb=false)
    dograb && grab(DATA_PATH)
    df = list_spectra()
    sort!(df, cols=:id, rev=true)

    if maximum(idx) < 9999
        # Treat as DataFrame Index
        data = load_spectra(df[:id][idx])
    else
        # Treat as Spectrum ID
        data = load_spectra(idx)
    end

    if length(idx) == 1
        # Make data iterable
        data = [data]
    end

    data = average.(data)

    # Empty Plot
    plot()
    for i = 1:length(data)
        if pixel
            x_binning = get_attribute(data[i], "x_binning")
            npixel = 512 / x_binning
            x = 1:npixel
            xlabelstr = "Pixel Number"
        elseif sfwl
            x = get_wavelength(data[i])
            xlabelstr = "Wavelength in nm"
        else
            x = get_ir_wavelength(data[i])
            xlabelstr = "IR Wavelength in nm"
        end

        name = get_attribute(data[i], "name")
        plot!(x, data[i].s, label=name)
        xlabel!(xlabelstr)
        ylabel!("Counts per Second")
    end
    gui()
    return nothing
end

"""
Process pump probe spectra. The data is splitted in the middle of the camera.

Keyword Arguments:

* `eventwidth`: Max. width of events that are being removed (default 8)
* `savepath`:   Path were the processed data is stored (default is the folder where the raw data is loaded from)
* `wlrange`:    Figures are cropped to these xlims (defaults to full range of available values) eg.: `wlrange=(3300, 3600)`
"""
function pump_probe_process(idarray;
    eventwidth=8,
    savepath="default",
    wlrange="default")

    print("Loading Data...")
    data = SFGTools.load_spectra(idarray)
    SFGTools.rm_events!.(data)
    data = SFGTools.average.(data)

    if savepath == "default"
        savepath = splitdir(splitdir(SFGTools.getdir(data[1].id))[1])[1]
        savepath = joinpath(savepath, "processed")
    end
    !isdir(savepath) && mkdir(savepath)

    dltime = SFGTools.get_pump_delay.(data)
    λ = SFGTools.get_ir_wavelength(data[1])
    ν = SFGTools.get_ir_wavenumber(data[1])
    name = SFGTools.get_attribute(data[1], "name")

    # Sort with respect to delay times
    p = sortperm(dltime)
    dltime = dltime[p]
    data = data[p]
    println("\tDone.")

    # Check if all delays are unique
    length(unique(dltime)) != length(dltime) && warn("Not all delay times are unique!")

    # Get Wavelength range
    wlrange == "default" && (wlrange = (minimum(λ), maximum(λ)))

    # Make seperate spectra for upper and lower part of CCD camera for all time
    # delays.
    print("Generating Spectra...")
    npixels = size(data[1].s)
    spectra = zeros(npixels[1],2)
    split = npixels[2] ÷ 2
    for d in data
        # Upper part
        spectra[:,1] += squeeze(sum(d.s[:,1:split], 2), 2)
        # Lower part
        spectra[:,2] += squeeze(sum(d.s[:,split+1:end], 2), 2)
    end
    # Normalize
    spectra /= length(data)

    plotarray = Plots.Plot[]
    p1 = Plots.plot(λ, spectra, labels=["A", "B"])
    Plots.title!(name * "/nSpectra Total")
    Plots.ylabel!("Counts [1/s]")
    Plots.xlabel!("IR Wavelength [nm]")
    Plots.xlims!(wlrange)
    push!(plotarray, p1)
    Plots.savefig(joinpath(savepath, "spectra.png"))
    println("\tDone.")

    # Make Heatmap
    print("Generating Heatmaps...")
    hmap = zeros(length(data), npixels[1], 2)
    for i in eachindex(data)
        hmap[i,:,1] = data[i].s[:,1:split]
        hmap[i,:,2] = data[i].s[:,split+1:end]
    end
    p2 = Plots.heatmap(λ, dltime, hmap[:,:,1], clim=(minimum(hmap), maximum(hmap)))
    Plots.title!(name * "\nA")
    Plots.xlabel!("IR Wavelength [nm]")
    Plots.ylabel!("Pump Delay [ps]")
    Plots.xlims!(wlrange)
    push!(plotarray, p2)
    Plots.savefig(joinpath(savepath, "heatmap_A.png"))

    p3 = heatmap(λ, dltime, hmap[:,:,2], clim=(minimum(hmap), maximum(hmap)))
    Plots.title!(name * "\nB")
    Plots.xlabel!("IR Wavelength [nm]")
    Plots.ylabel!("Pump Delay [ps]")
    Plots.xlims!(wlrange)
    push!(plotarray, p3)
    Plots.savefig(joinpath(savepath, "heatmap_B.png"))
    println("\tDone.")

    # Animated Spectra
    print("Animating Spectra")
    labels=("A", "B")
    anim = Plots.@animate for i = 1:size(hmap, 1)
        print(".")

        p1 = Plots.plot(λ, hmap[i,:,:], labels=labels, ylim=(minimum(hmap), maximum(hmap)))
        Plots.title!("Delay = $(round(dltime[i], 1)) ps")
        Plots.ylabel!("Counts [1/s]")
        Plots.xlabel!("IR Wavelength [nm]")
        Plots.xlims!(wlrange)

        p2 = Plots.heatmap(λ, dltime, hmap[:,:,1], clim=(minimum(hmap), maximum(hmap)))
        Plots.title!(name * "\nA")
        Plots.xlabel!("IR Wavelength [nm]")
        Plots.ylabel!("Pump Delay [ps]")
        Plots.scatter!([wlrange[1], wlrange[end]], [dltime[i], dltime[i]], markersize=10, legend=false, color=:white)
        Plots.xlims!(wlrange)

        p3 = heatmap(λ, dltime, hmap[:,:,2], clim=(minimum(hmap), maximum(hmap)))
        Plots.title!(name * "\nB")
        Plots.xlabel!("IR Wavelength [nm]")
        Plots.ylabel!("Pump Delay [ps]")
        Plots.scatter!([wlrange[1], wlrange[end]], [dltime[i], dltime[i]], markersize=10, legend=false, color=:white)
        Plots.xlims!(wlrange)

        Plots.plot(p1, p2, p3, layout=(3,1), size=(600,1200))
    end
    Plots.gif(anim, joinpath(savepath, "spectra.gif"), fps=2)
    println("\tDone.")

    return spectra, plotarray
end

function cal_spectrometer(idarray)

    data = SFGTools.load_spectra(idarray)
    data = SFGTools.average.(data)
    names = get_attribute(data, "name")
    splitted_names = split.(names, "_")

    λcal = zeros(length(data))
    locs = fill("", length(data))
    for i in eachindex(splitted_names)
        λcal[i] = splitted_names[i][end-1] |> parse
        locs[i] = splitted_names[i][end]
    end

    # Check if the data is suitable
    for i = 1:length(data)
        length(data[i].s) != 512 && error("All spectra need to have 512 data points")
        !any(["center", "right", "left"] .== locs[i]) && error("Filenames have to have the form \"*_334.15_left\" etc. (see help)")
    end

    # Remove everything but the peaks
    for i in eachindex(data)
        if locs[i] == "center"
            data[i].s[[1:236..., 276:end...]] = NaN
        elseif locs[i] == "left"
            data[i].s[20:end] = NaN
        elseif locs[i] == "right"
            data[i].s[1:end-20] = NaN
        end
    end

    peakpositions = fill(0, length(data))
    for i in 1:length(data)
        peakpositions[i] = findmax(data[i].s)[2]
    end
    peakpositions = reshape(peakpositions, (3, length(data) ÷ 3))'
    peakpositions = hcat(peakpositions[:,2], peakpositions[:,1], peakpositions[:,3])

    λ_spectrometer = Float64[]
    for (i, d) in enumerate(data)
        λ = get_attribute(d, "spectrometer_wavelength")
        push!(λ_spectrometer, λ)
    end
    λ_spectrometer = reshape(λ_spectrometer, (3, 7))'
    λ_spectrometer = hcat(λ_spectrometer[:,2], λ_spectrometer[:,1], λ_spectrometer[:,3])

    Δp = abs.(diff(peakpositions, 2))
    Δλ = abs.(diff(λ_spectrometer, 2))
    nmperpixel = Δλ ./ Δp

    mean_nmperpixel = mean(nmperpixel, 2)
    mean_nmperpixel = squeeze(mean_nmperpixel, 2)

    model(x, p) = p[1] * x.^2 + p[2] * x + p[3]
    fit = curve_fit(model, λ_spectrometer[:,2], mean_nmperpixel, [-3e-8, 5e-6, 0.02])
    xfit = linspace(330, 600, 100)
    yfit = model(xfit, fit.param)

    p1 = scatter(λ_spectrometer[:,2], mean_nmperpixel)
    plot!(xfit, yfit)
    xlabel!("Central Spectrometer Wavelength in nm")
    ylabel!("Wavelength Density [nm/pixel]")

    # the real center is in between pixels 256 and 257!
    Δp0 = peakpositions[:,2] - 256.5

    offset_fit = curve_fit(model, λ_spectrometer[:,2], Δp0, [-3e-4, 0.25, -32.0])
    offset_yfit = model(xfit, offset_fit.param)

    p2 = scatter(λ_spectrometer[:,2], Δp0)
    plot!(xfit, offset_yfit)
    xlabel!("Central Spectrometer Wavelength in nm")
    ylabel!("Pixel Offset")

    plot(p1, p2)
    gui()

    info("Add the following code to the `get_wavelength` function at the appropriate location in the SFGTools package:")

    println("
    # Calibration parameters were determined on $(now())
    # Pixel offset
    Δp = $(offset_fit.param[1]) * λ0^2 + $(offset_fit.param[2]) * λ0 - $(offset_fit.param[3])
    # Wavelenths per pixel
    pλ = $(fit.param[1]) * λ0^2 + $(fit.param[2]) * λ0 + $(fit.param[3])
    ")

    info("The right location is inside the if-elseif-else statement. Change the last else to:")

    println("elseif date < DateTime(\"$(now())\")")

    info("and add a new final else statement with the above code that contains the calibration parameters.")
    return nothing
end

end #module

max_stars_repo_name: PhysChemUniDue/DataScripts

max_issues_repo_name: PhysChemUniDue/DataScripts

max_forks_repo_name: PhysChemUniDue/DataScripts


text:
export gather!

import MPI

@doc """
    gather!(A, A_global)
    gather!(A, A_global; root=0)

Gather a CPU-array `A` from each member of the Cartesian grid of MPI processes into a one large CPU-array `A_global` on the root process (default: `0`).

!!! note "Memory usage note"
    `gather!` allocates at first call an internal buffer of the size of `A_global` and keeps it alive until `finalize_global_grid` is called. A (re-)allocation occurs only if `gather!` is called with a larger `A_global` than in any previous call since the call to `init_global_grid`. This is an optimisation to minimize (re-)allocation, which is very important as `gather!` is typically called in the main loop of a simulation and its performance is critical for the overall application performance.
"""
gather!

let
    global gather!, free_gather_buffer
    A_all_buf = zeros(0);

    "Free the buffer used by gather!."
    function free_gather_buffer()
        A_all_buf = nothing;
        GC.gc();
        A_all_buf = zeros(0);
    end

    function gather!(A::Array{T}, A_global::Union{Array{T}, Nothing}; root::Integer=0) where T <: GGNumber
        check_initialized();
        cart_gather!(A, A_global, me(), global_grid().dims, comm(); root=root);
        return nothing
    end

    function cart_gather!(A::Array{T}, A_global::Union{Array{T}, Nothing}, me::Integer, dims::Array{T2}, comm_cart::MPI.Comm; root::Integer=0, tag::Integer=0, ndims::Integer=NDIMS_MPI) where T <: GGNumber where T2 <: Integer
        nprocs = prod(dims);
        if me != root
            req = MPI.REQUEST_NULL;
            req = MPI.Isend(A, root, tag, comm_cart);
            MPI.Wait!(req);
        else # (me == root)
            A_global === nothing && error("The input argument A_global can't be `nothing` on the root")
            if length(A_global) != nprocs*length(A) error("The input argument A_global must be of length nprocs*length(A)") end
            if (eltype(A_all_buf) != T)
                A_all_buf = reinterpret(T, A_all_buf);
            end
            if length(A_all_buf) < nprocs*length(A)                       # Allocate only if the buffer is not large enough
                free_gather_buffer();                                     # Free the memory of the old buffer immediately as it can typically go up to the order of the total available memory.
                A_all_buf = zeros(T, Int(ceil(nprocs*length(A)/GG_ALLOC_GRANULARITY))*GG_ALLOC_GRANULARITY);  # Ensure that the amount of allocated memory is a multiple of GG_ALLOC_GRANULARITY*sizeof(T). So, we can always correctly reinterpret A_all_buf even if next time sizeof(T) is greater.
            end
            A_all_flat = view(A_all_buf,1:nprocs*length(A));              # Create a 1D-view on the amount of memory needed from A_all_buf.
            reqs = fill(MPI.REQUEST_NULL, nprocs);
            for p in [0:root-1; root+1:nprocs-1]
                cs = Cint[-1,-1,-1];
                MPI.Cart_coords!(comm_cart, p, cs);
                offset = cs[1]*length(A) + cs[2]*dims[1]*length(A) + cs[3]*dims[1]*dims[2]*length(A)
                A_c = view(A_all_flat, 1+offset:length(A)+offset);
                reqs[p+1] = MPI.Irecv!(A_c, p, tag, comm_cart);           # Irev! requires a contigous (SubArray) buffer (that is not both reshaped and reinterpreted)...
            end
            cs = MPI.Cart_coords(comm_cart);
            A_all = reshape(A_all_flat, (length(A), dims[1], dims[2], dims[3])); # Create a 4D-view on the amount of memory needed from A_all_buf.
            A_all[:,cs[1]+1,cs[2]+1,cs[3]+1] .= A[:];
            if (nprocs>1) MPI.Waitall!(reqs); end
            nx, ny, nz = size(view(A,:,:,:));
            for cz = 0:size(A_all,4)-1, cy = 0:size(A_all,3)-1, cx = 0:size(A_all,2)-1
                A_global[(1:nx).+cx*nx, (1:ny).+cy*ny, (1:nz).+cz*nz] .= reshape(A_all[:,cx+1,cy+1,cz+1],nx,ny,nz); # Store the data at the right place in A_global (works for 1D-3D, e.g. if 2D: nz=1, cz=0...)
            end
        end
    end
end

ext: jl


text:
import StaticArrays: SVector, @SVector

function distsq(A::SVector{2, T}, B::SVector{2, T}) where T <: Real
    return (A[1]-B[1])^2 + (A[2]-B[2])^2
end

function grid_idxs(point:: SVector{2, T}, cellsize::T) where T <: Real
    """ Returns cell indices, (i,j)=(row, col) for  `point` """
    return ceil(Int, point[2] / cellsize) + 1, ceil(Int, point[1] / cellsize) + 1
end

function ispoint_inrectangle(point::SVector{2, T}, rlims::Tuple{T, T}) where T <: Real
    """
    Check whether a point is within the rectangle
    [0, rlims[1]] x [0, rlims[2]];
    here rlims =(width, height)
    """
    return (0 <= point[1] <= rlims[1]) && (0 <= point[2] <= rlims[2])
end

function unif_sampling_annulus(center::SVector{2, T}, inner_rad::T, outer_rad::T;
                               k=30) where T <: Real
    """
    Samples  uniformly k points from the annulus of center, `center`,
    and radii, inner_radius, outer_radius
    """
    rad = sqrt.(inner_rad^2  .+ (outer_rad^2 - inner_rad^2)*rand(k))
    angle = 2*pi * rand(k)
    x, y = center[1] .+ rad .* cos.(angle), center[2] .+ rad .*sin.(angle)
    return [@SVector([xi, yi]) for (xi, yi) in zip(x,y)]
end

function Poisson_disk_sampling(;width =1.0, height=1.0, radius=0.05, k=30)
    """
    Bridson Algorithm to generate  Poisson disk samples
    within a 2D rectangle, [0, width] x [0, height]
    Ref: R Bridson, Fast Poisson disk sampling in arbitrary dimensions,
    in  SIGGRAPH '07: ACM SIGGRAPH 2007 sketches
    https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf
    """
    function  is_accepted(point:: SVector{2, T})  where T <: Real
        """ Check whether the point, `point` can be added to the active_list"""

        i, j = grid_idxs(point, cellsize)
        m, n = maximum([i-2, 1]), maximum([j-2, 1])
        p, q = minimum([i+2, grid_rows]), minimum([j+2, grid_cols])
        for k in n:q
            for l in m:p
                gpoint = grid[(k-1) * grows + l]
                if ismissing(gpoint)
                    continue
                end
                if distsq(point, gpoint) <= r2
                    return false
                end
            end
        end
        return true
    end

    cellsize = radius/sqrt(2)
    grid_cols = ceil(Int, width / cellsize) + 1
    grid_rows = ceil(Int, height / cellsize) + 1

    grid = Vector{Any}(missing, grid_rows*grid_cols)
    r2 = radius*radius
    active_list = SVector{2, Float64}[]

    #generate a random point within the rectangle of interest
    # and  insert it into the active_list and grid
    point = @SVector([width*rand(1)[1], height*rand(1)[1]])
    push!(active_list, point)
    i, j = grid_idxs(point, cellsize)
    grid[(j-1)*grid_rows+i] = point

    while length(active_list) > 0
        i = rand(1:length(active_list), 1)[1]
        pt = active_list[i]
        active_list  =  active_list[1:end .!= i] #removes the point in pos i
        rpoints = unif_sampling_annulus(pt, radius, 2*radius; k=k)
        for rp in rpoints
            if !ispoint_inrectangle(rp, (width, height)) || !is_accepted(rp)
                continue
            end
            push!(active_list, rp)
            i, j = grid_idxs(rp, cellsize)
            grid[(j-1)*grid_rows+i] = rp
        end
    end
    return filter(x->!ismissing(x), grid);
end

max_stars_repo_name: empet/PoissonDisk

max_issues_repo_name: empet/PoissonDisk

max_forks_repo_name: empet/PoissonDisk


text:
# Tests for the Diffusion module
#

function test_diffusion_basic()
    model = Diffusion.Model(N=4, L=2, K=0.1)
    model.parameters.K == 0.1
end

function test_diffusion_set_c()
    model = Diffusion.Model(N=4, L=2, K=0.1)
    c0 = 1:4
    model.solution.c = c0
    model.solution.c.data[1:model.grid.N] == c0
end

function test_diffusive_flux(stepper=:ForwardEuler; top_flux=0.3, bottom_flux=0.13, N=10)
    model = Diffusion.Model(N=N, L=1, K=1.0, stepper=stepper)
    model.bcs.c.top = FluxBoundaryCondition(top_flux)
    model.bcs.c.bottom = FluxBoundaryCondition(bottom_flux)

    C₀ = integral(model.solution.c)
    C(t) = C₀ - (top_flux - bottom_flux) * t

    dt = 1e-6
    iterate!(model, dt, 10)

    return C(time(model)) ≈ integral(model.solution.c)
end

function test_diffusion_cosine(stepper=:ForwardEuler)
    model = Diffusion.Model(N=100, L=π/2, K=1.0, stepper=stepper)
    z = model.grid.zc

    c_init(z) = cos(2z)
    c_ans(z, t) = exp(-4t) * c_init(z)

    model.solution.c = c_init

    dt = 1e-3
    iterate!(model, dt)

    # The error tolerance is a bit arbitrary.
    norm(c_ans.(z, time(model)) .- data(model.solution.c)) < model.grid.N*1e-6
end

function test_diffusion_cosine_run_until(stepper=:ForwardEuler)
    model = Diffusion.Model(N=100, L=π/2, K=1.0, stepper=stepper)
    z = model.grid.zc

    c_init(z) = cos(2z)
    c_ans(z, t) = exp(-4t) * c_init(z)

    model.solution.c = c_init

    dt = 1e-3
    tfinal = 3dt/2
    run_until!(model, dt, tfinal)

    # The error tolerance is a bit arbitrary.
    norm(c_ans.(z, tfinal) .- data(model.solution.c)) < model.grid.N*1e-6
end

function test_advection(stepper=:ForwardEuler)
    L = 1.0
    W = -1.0
    δ = L/10
    h = L/2

    model = Diffusion.Model(N=100, L=L, K=0.0, W=W, stepper=stepper,
                            bcs = Diffusion.BoundaryConditions(FieldBoundaryConditions(
                                    GradientBoundaryCondition(0.0),
                                    GradientBoundaryCondition(0.0))))

    c_gauss(z, t) = exp( -(z-W*t)^2 / (2*δ^2) )
    c₀(z) = c_gauss(z+h, 0)
    model.solution.c = c₀

    iterate!(model, 1e-3, 10)

    c_current(z) = c_gauss(z+h, time(model))
    c_answer = CellField(model.grid)
    set!(c_answer, c_current)

    norm(data(c_answer) .- data(model.solution.c)) < 0.05
end

function test_damping(stepper=:ForwardEuler)
    L = 1.0
    μ = 0.1

    model = Diffusion.Model(N=3, L=L, K=0.0, μ=μ, stepper=stepper,
                            bcs = Diffusion.BoundaryConditions(FieldBoundaryConditions(
                                    GradientBoundaryCondition(0.0),
                                    GradientBoundaryCondition(0.0))))

    c_damp(t) = exp(-μ*t)
    c₀(z) = c_damp(0)
    model.solution.c = c₀

    iterate!(model, 1e-3, 100)

    c_current(t) = c_damp(time(model))
    c_answer = CellField(model.grid)
    set!(c_answer, c_current)

    all(abs.(data(c_answer) .- data(model.solution.c)) .< 1e-6)
end

@testset "Diffusion" begin
    @test test_diffusion_basic()
    @test test_diffusion_set_c()

    for stepper in steppers
        @test test_diffusion_cosine(stepper)
        @test test_diffusion_cosine_run_until(stepper)
        @test test_diffusive_flux(stepper, top_flux=0, bottom_flux=0)
        @test test_diffusive_flux(stepper)
        @test test_advection(stepper)
        @test test_damping(stepper)
    end
end


text:
##################################################
# Hamiltonians
##################################################

"""
    overlap([T::Type], f, args...; kwargs...)

Construct the matrix corresponding to the overlap operator of an atom,
bond or molecule.  `f` is a function that accepts the atom types and
distances between the atoms to compute the overlap matrix element.

"""
overlap(f, args...; kwargs...) = overlap(Complex{Float64}, f, args...; kwargs...)

function overlap(T::Type, f, atom::Atom; kwargs...)
    N = countorbitals(atom)
    H = fill(zero(T), N, N)
    U = localbasis(atom)
    for (i, s1) in enumerate(atom.orbitals)
        H[i,i] += f(chemical_symbol(atom) => shortform(s1))
    end
    return U'H*U
end

function overlap(T::Type, f, atom1::Atom, atom2::Atom; fixed_distance=nothing, kwargs...)
    N1 = countorbitals(atom1)
    N2 = countorbitals(atom2)
    U1 = localbasis(atom1)
    U2 = localbasis(atom2)
    H = fill(zero(T), N1, N2)

    if fixed_distance === nothing
        r = atom1.position - atom2.position
    else
        r = fixed_distance .* normalize(atom1.position - atom2.position)
    end

    for (j, s2) in enumerate(atom2.orbitals), (i, s1) in enumerate(atom1.orbitals)
        H[i,j] += f(r, chemical_symbol(atom1) => shortform(s1), chemical_symbol(atom2) => shortform(s2))
    end
    return U1'H*U2
end

function overlap(T::Type, f, mol::Molecule; kwargs...)
    N = countorbitals(mol)
    H = fill(zero(T), N, N)

    i = 1
    for atom1 in mol.atoms
        j = 1
        Ni = countorbitals(atom1)
        for atom2 in mol.atoms
            Nj = countorbitals(atom2)
            if atom1 === atom2
                H[i:i+Ni-1, j:j+Nj-1] += overlap(T, f, atom1; kwargs...)
            elseif Bond(atom1, atom2) in mol.bonds
                H[i:i+Ni-1, j:j+Nj-1] += overlap(T, f, atom1, atom2; kwargs...)
            end
            j += Nj
        end
        i += Ni
    end
    return H
end

function overlap(T::Type, f, mol1::Molecule, mol2::Molecule, bonds::Set{<:Bond}; kwargs...)
    N1 = countorbitals(mol1)
    N2 = countorbitals(mol2)
    H = fill(zero(T), N1, N2)

    i = 1
    for atom1 in mol1.atoms
        j = 1
        Ni = countorbitals(atom1)
        for atom2 in mol2.atoms
            Nj = countorbitals(atom2)
            if Bond(atom1, atom2) in bonds
                H[i:i+Ni-1, j:j+Nj-1] += overlap(T, f, atom1, atom2; kwargs...)
            end
            j += Nj
        end
        i += Ni
    end
    return H
end


text:
"""
    Circuit(n_nodes::Int)

A Hamiltonian circuit. If the vector `x` is constrained within a `Circuit(n)`,
each `x[i]` denotes the next node in the graph, for `i ∈ [1, n]`.

The considered graph is an undirected complete graph with `n` nodes.

Also called `cycle` or `atour`.
"""
struct Circuit <: MOI.AbstractVectorSet
    n_nodes::Int
end

MOI.dimension(set::Circuit) = set.n_nodes

"""
    CircuitPath(n_nodes::Int)

A Hamiltonian circuit. If the vectors `x` and `y` are constrained within a
`CircuitPath(n)`, each `x[i]` denotes the next node in the graph, for
`i ∈ [1, n]`. The last `n` variables denote the order in which the nodes are
visited, i.e. `y[1]` is the first visited node (1 by convention), `y[2]` is
the next node in the path, etc.

The considered graph is an undirected complete graph with `n` nodes.
"""
struct CircuitPath <: MOI.AbstractVectorSet
    n_nodes::Int
end

MOI.dimension(set::CircuitPath) = 2 * set.n_nodes

"""
    WeightedCircuit{T <: Real}(n_nodes::Int, cost_matrix::AbstractMatrix{T})

A Hamiltonian circuit. If the vector `x` and the scalar `c` are constrained
within a `WeightedCircuit(n, cost_matrix)`, each `x[i]` denotes the next node
in the graph, for `i ∈ [1, n]`. `c` is the total cost of the circuit, defined as:

``c = \\sum_{i=1}^n \\mathtt{cost\\_matrix}_{i, x[i]}``

The considered graph is an undirected complete graph with `n` nodes.
"""
struct WeightedCircuit{T <: Real} <: MOI.AbstractVectorSet
    n_nodes::Int
    cost_matrix::AbstractMatrix{T}
end

MOI.dimension(set::WeightedCircuit{T}) where {T} = set.n_nodes + 1
function copy(set::WeightedCircuit{T}) where {T}
    return WeightedCircuit(set.n_nodes, copy(set.cost_matrix))
end
function Base.:(==)(x::WeightedCircuit{T}, y::WeightedCircuit{T}) where {T}
    return x.n_nodes == y.n_nodes && x.cost_matrix == y.cost_matrix
end

"""
    WeightedCircuitPath(n_nodes::Int, cost_matrix::AbstractMatrix{T})

A Hamiltonian circuit. If the vectors `x` and `y` and the scalar `c` are
constrained within a `CircuitPath(n)`, each `x[i]` denotes the next node in the graph, for
`i ∈ [1, n]`. The next `n` variables denote the order in which the nodes are
visited, i.e. `y[1]` is the first visited node (1 by convention), `y[2]` is
the next node in the path, etc. `c` is the total cost of the circuit, defined as:

``c = \\sum_{i=1}^n \\mathtt{cost\\_matrix}_{i, x[i]}``

The considered graph is an undirected complete graph with `n` nodes.
"""
struct WeightedCircuitPath{T <: Real} <: MOI.AbstractVectorSet
    n_nodes::Int
    cost_matrix::AbstractMatrix{T}
end

MOI.dimension(set::WeightedCircuitPath{T}) where {T} = 2 * set.n_nodes + 1
function copy(set::WeightedCircuitPath{T}) where {T}
    return WeightedCircuitPath(set.n_nodes, copy(set.cost_matrix))
end
function Base.:(==)(
    x::WeightedCircuitPath{T},
    y::WeightedCircuitPath{T},
) where {T}
    return x.n_nodes == y.n_nodes && x.cost_matrix == y.cost_matrix
end

# isbits types, nothing to copy
function copy(set::Union{Circuit, CircuitPath})
    return set
end


text:
# ---
# title: 698. Partition to K Equal Sum Subsets
# id: problem698
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Recursion
# link: <https://leetcode.com/problems/partition-to-k-equal-sum-subsets/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and a positive integer `k`, find whether
# it's possible to divide this array into `k` non-empty subsets whose sums are
# all equal.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [4, 3, 2, 3, 5, 2, 1], k = 4
#     Output: True
#     Explanation: It's possible to divide it into 4 subsets (5), (1, 4), (2,3), (2,3) with equal sums.
#
#
#
#
# **Note:**
#
#   * `1 <= k <= len(nums) <= 16`.
#   * `0 < nums[i] < 10000`.
#
#
## @lc code=start
using LeetCode

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


text:
module Jl2Cpp
    DLL_EXPORT = "#ifdef _MSC_VER\n__declspec(dllexport)\n#endif\nextern \"C\""
    #CC = "g++"
    CC = haskey(ENV, "JL2CPP_CC")? ENV["JL2CPP_CC"]: "D:\\LLVM\\bin\\clang++.exe"
    #C_OPT = ["-std=c++14", """-I$(joinpath(Base.source_dir(), "includes"))""", "-shared", "-Ofast", "-g"]
    C_OPT = haskey(ENV, "JL2CPP_C_OPT")? [split(ENV["JL2CPP_C_OPT"])..., """-I$(joinpath(Base.source_dir(), "includes"))"""]: ["-std=c++14", """-I$(joinpath(Base.source_dir(), "includes"))""", "-shared", "-Ofast", "-g"]
    CPP_PREDEFS = []
    CPP_INCLUDES = ["#include<cstdint>", "#include<tuple>", "#include<algorithm>", "#include<complex>", "#include<cstdio>", "#include<cfenv>", "#include<cassert>", "#include<jl2cpp.hpp>", "#include<Eigen/Dense>"]
    CPP_USINGS = ["using namespace Eigen;"]

    """
    Type to pass matrixes or arrays to C++ function
    """
    mutable struct CMatrix{T}
        ptr::Ptr{T}
        rsize::Csize_t
        csize::Csize_t
    end
    export CMatrix
    mutable struct CTuple
        ptr::Ptr{Void}
        size::Csize_t
    end
    export CTuple

    function ConvertToArray{T}(m::CMatrix{T})
        # ret = (m.csize == 1)? Array{T, 1}(m.rsize): Array{T, 2}(m.rsize, m.csize)
        unsafe_wrap((m.csize == 1)? Array{T, 1}: Array{T, 2}, m.ptr, (m.csize == 1)?m.rsize:(m.rsize, m.csize), false)
    end
    export ConvertToArray

    function ConvertToTuple(ty::DataType, from::CTuple)
        assert(ty.name.name == :Tuple)
        ret = :((1,))
        ret.args = map(x->:(zero($x)), ty.parameters)
        ret = eval(ret)
        offset = 0
        for i=1:length(ret)
            unsafe_copy!(
                convert(Ptr{fieldtype(ty, i)}, pointer_from_objref(ret) + fieldoffset(ty, i)),
                convert(Ptr{fieldtype(ty, i)}, from.ptr + offset),
                1
            )
            offset += sizeof(fieldtype(ty, i))
            if offset > from.size; error("wrong tuple size"); end
        end
        ret
    end
    export ConvertToTuple

    mutable struct GenCppCodeInfo
        indent :: Int64
        userDefFunc::Set{String}
        argArrays::Set{Symbol}  # arrays of user-defined function are safe to get pointer, so use them to speedup
    end
    GenCppCodeInfo() = GenCppCodeInfo(0, Set([]), Set([]))

    function putInfo(str)
        if !haskey(ENV, "JL2CPP_QUIET")
            println(str)
        end
    end
    function IsType(x)
        if x isa DataType; return true; end
        if x isa Expr && x.head == :call && x.args[1] in [:eltype]; return true; end # for functions that return type
        try
            x isa Symbol && eval(x) isa DataType
        catch
            false
        end
    end
    function GetType(x, info)
        if x isa DataType; return ConvertType(x); end
        if x isa Expr && x.head == :call && x.args[1] in [:eltype]; return GenerateCppCode(x, info); end # for functions that return type
        if x isa Symbol; return ConvertType(eval(x)); end
    end
    RECIPES = []
    function __init__()
        for idx in find(x->x=="-r", ARGS)
            rn = ARGS[idx + 1]
            putInfo("recipe: $rn")
            push!(RECIPES, evalfile(joinpath(@Base.__DIR__, "recipe", rn * ".jl")))
        end
    end

    function GenerateWrapperFunction(ex::Expr, info)::String
        if ex.args[1].head != :(::)
            error("functions must be typed")
        end
        name = ex.args[1].args[1].args[1]
        args = ex.args[1].args[1].args[2:end]
        retty = ex.args[1].args[2] |> eval
        rettycpp = if retty.name.name == :Array
            "jl2cpp::matrix"
        elseif retty.name.name == :Tuple
            "jl2cpp::tuple"
        else
            ConvertType(retty)
        end
        argstrs = String[]
        initstrs = String[]
        argnames = String[]
        info.argArrays = Set([])
        push!(initstrs, "initParallel();")
        for a = args
            if typeof(a) != Expr || a.head != :(::)
                error("function arguments must be typed")
            end
            cname = a.args[1]
            push!(argnames, cname)
            ty = eval(a.args[2])
            if ty.name.name == :Array
                cty = ty.parameters[1] |> ConvertType
                if ty.parameters[2] == 1
                    push!(argstrs, "jl2cpp::matrix *$(cname)_jl2cpp_arg")
                    push!(initstrs, repeat("\t", info.indent + 1) * "Map<Matrix<$cty, Dynamic, 1> > $cname(($cty *)$(cname)_jl2cpp_arg->ptr, $(cname)_jl2cpp_arg->rsize);")
                    #push!(initstrs, repeat("\t", info.indent + 1) * "$cty *$(cname)_jl2cpp_ptr = ($cty *)$(cname)_jl2cpp_arg->ptr;")
                    #push!(initstrs, repeat("\t", info.indent + 1) * "printf(\"size: %d\\n\", $(cname)_jl2cpp_arg->rsize);")
                elseif ty.parameters[2] == 2
                    push!(argstrs, "jl2cpp::matrix *$(cname)_jl2cpp_arg")
                    push!(initstrs, repeat("\t", info.indent + 1) * "Map<Matrix<$cty, Dynamic, Dynamic> > $cname(($cty *)$(cname)_jl2cpp_arg->ptr, $(cname)_jl2cpp_arg->rsize, $(cname)_jl2cpp_arg->csize);")
                    #push!(initstrs, repeat("\t", info.indent + 1) * "$cty *$(cname)_jl2cpp_ptr = ($cty *)$(cname)_jl2cpp_arg->ptr;")
                else
                    error("3+ dimension array is currently not supported.")
                end
            else
                cty = ty |> ConvertType
                push!(argstrs, "$cty $cname")
            end
        end
        argstr = join(argstrs, ", ")
        initstr = join(initstrs, "\n")
        retstr = ""
        if retty.name.name == :Void
            retstr = "$name($(join(argnames, ", ")));"
        elseif retty.name.name == :Tuple
            retstr = "auto jl2cpp_ret_tmp = $name($(join(argnames, ", ")));\n"
            retstr *= repeat("\t", info.indent + 1) * "void* jl2cpp_ret_ptr = malloc(jl2cpp::SumTupleSize<decltype(jl2cpp_ret_tmp)>::value);\n"
            retstr *= repeat("\t", info.indent + 1) * "std::size_t jl2cpp_ret_offset = 0;\n"
            i = 0
            for t in retty.parameters
                tcppty = ConvertType(t)
                retstr *= repeat("\t", info.indent + 1) * "*($tcppty *)((std::uint8_t *)jl2cpp_ret_ptr + jl2cpp_ret_offset) = std::get<$i>(jl2cpp_ret_tmp);\n"
                retstr *= repeat("\t", info.indent + 1) * "jl2cpp_ret_offset += sizeof($tcppty);\n"
                i += 1
            end
            retstr *= repeat("\t", info.indent + 1) * "return jl2cpp::tuple{jl2cpp_ret_ptr, jl2cpp_ret_offset};\n"
        else
            retstr = "return jl2cpp::calc_retval($name($(join(argnames, ", "))));"
        end

        body = "{\n" * initstr * (isempty(initstrs)? "": "\n") * repeat("\t", info.indent + 1) * retstr * "\n}"
        "$DLL_EXPORT $rettycpp $(name)_jl2cpp_wrapper($argstr)\n$body\n"
    end

    """
    Convert Julia type into C++ type.
    """
    function ConvertType(ty :: DataType)
        sym = ty.name.name
        if sym == :Int64
            "int64_t"
        elseif sym == :Int32
            "int32_t"
        elseif sym == :UInt64
            "uint64_t"
        elseif sym == :UInt32
            "uint32_t"
        elseif sym == :String
            "std::string"
        elseif sym == :Float64
            "double"
        elseif sym == :Float32
            "float"
        elseif sym == :Complex64
            "std::complex<float>"
        elseif sym == :Complex128
            "std::complex<double>"
        elseif sym == :Complex
            t = ConvertType(ty.parameters[1])
            "std::complex<$t>"
        elseif sym == :Void
            "void"
        elseif sym == :Tuple
            "std::tuple<"* join(map(x->if x.name.name in [:Array, :Tuple]; error("nested tuple is not supported");else; ConvertType(x);end, ty.parameters), ", ") * ">"
        elseif sym == :Array
            if ty.parameters[1] in [:Tuple, :Array]; error("nested array is not supported");end
            t = ConvertType(ty.parameters[1])
            if ty.parameters[2] == 1
                "Matrix<$t, Dynamic, 1>"
            elseif ty.parameters[2] == 2
                "Matrix<$t, Dynamic, Dynamic>"
            else
                error("3+ dimension array is not supported.")
            end
        elseif sym == :Char
            "char"
        else
            error("unknown type: $sym")
        end
    end

    """
    Convert Julia function into C++ form
    """
    function ConvertFunction(name :: String, args, info) :: String
        if name == "length"
            "$(GenerateCppCode(args[1], info)).size()"
        elseif name == "size"
            if length(args) != 2
                "std::forward_as_tuple($(GenerateCppCode(args[1], info)).rows(), $(GenerateCppCode(args[1], info)).cols())"
            elseif contains(GenerateCppCode(args[2], info), "1")
                "$(GenerateCppCode(args[1], info)).rows()"
            elseif contains(GenerateCppCode(args[2], info), "2")
                "$(GenerateCppCode(args[1], info)).cols()"
            else
                error("3+ dimension array is not supported.")
            end
        elseif name == "eye"
            "MatrixXd::Identity($(GenerateCppCode(args[1], info)), $(length(args) < 2? 1: GenerateCppCode(args[2], info))).eval()"
        elseif name == "zeros"
            if IsType(args[1])
                a = GetType(args[1], info)
                "Matrix<$a, Dynamic, Dynamic>::Zero($(GenerateCppCode(args[2], info)), $(length(args) < 3? 1: GenerateCppCode(args[3], info))).eval()"
            else
                a = GenerateCppCode(args[1], info)
                "MatrixXd::Zero($a, $(length(args) < 2? 1: GenerateCppCode(args[2], info))).eval()"
            end
        elseif name == "ones"
            if IsType(args[1])
                a = GetType(args[1], info)
                "Matrix<$a, Dynamic, Dynamic>::Ones($(GenerateCppCode(args[2], info)), $(length(args) < 3? 1: GenerateCppCode(args[3], info))).eval()"
            else
                a = GenerateCppCode(args[1], info)
                "MatrixXd::Ones($a, $(length(args) < 2? 1: GenerateCppCode(args[2], info))).eval()"
            end
        elseif name == "rand"
            #"(MatrixXd::Random($(GenerateCppCode(args[1], info)), $(GenerateCppCode(args[2], info))) + 1) / 2"
            # since arguments of `rand` is much complicated, the implementation is WIP
            if length(args) == 0
                "Matrix<double, 1, 1>::Random()[0]"
            elseif args[1] |>IsType
                typ = GetType(args[1], info)
                if args[1] in [:Float, :Float64, :Float32, :Float16, :Complex, :Complex128, :Complex64, :Complex32]
                    "((Matrix<$typ, 1, 1>::Random()[0] + ($typ)(1.0)) / 2)"
                else
                    "Matrix<$typ, 1, 1>::Random()[0]"
                end
            elseif isa(args[1], Expr) && args[1].head == :(:)
                error("not implemented yet")
            else
                # all arguments assumed dimension parameter
                if length(args) == 1
                    "jl2cpp::random($(GenerateCppCode(args[1], info)))"
                elseif length(args) == 2
                    "jl2cpp::random($(GenerateCppCode(args[1], info)),$(GenerateCppCode(args[2], info)))"
                else
                    error("3+ dimension is currently not supported")
                end
            end
        elseif name == "randn"
            "jl2cpp::randn($(join(map(x->GenerateCppCode(x, info), args), ", ")))" # TODO: implement type argument
        elseif name == "mean"
            "$(GenerateCppCode(args[1], info)).mean()"
        elseif name == "trace"
            "$(GenerateCppCode(args[1], info)).trace()"
        elseif name == "std"
            "jl2cpp::stdDev($(GenerateCppCode(args[1], info)))"
        elseif name == "sum"
            "$(GenerateCppCode(args[1], info)).sum()"
        elseif name == "eltype"
            "typename jl2cpp::ElType<decltype($(GenerateCppCode(args[1], info)))>::type::Scalar"
        elseif name == "zero"
            if IsType(args[1])
                a = GetType(args[1], info)
                "static_cast<$a>(0)"
            else
                a = GenerateCppCode(args[1], info)
                "static_cast<jl2cpp::ElType<decltype($a)>::type>(0)"
            end
        elseif name == "one"
            if IsType(args[1])
                a = GetType(args[1], info)
                "static_cast<$a>(1)"
            else
                a = GenerateCppCode(args[1], info)
                "static_cast<jl2cpp::ElType<decltype($a)>::type>(1)"
            end
        elseif name == "hex"
            "jl2cpp::hex($(GenerateCppCode(args[1], info)))"
        elseif name == "parse"
            if length(args) < 1 || !IsType(args[1])
                error("parsing julia expression is not supported")
            end
            if length(args) < 2
                error("too few argument in parse")
            end
            typ = args[1]|>eval|>ConvertType
            base = length(args) == 3? GenerateCppCode(args[3], info): 10
            "jl2cpp::parse<$typ>($(GenerateCppCode(args[2], info)), $base)"
        elseif name == "max"
            if length(args) == 2
                argstr = map(x->GenerateCppCode(x, info), args)
                "($(argstr[1]) > $(argstr[2])? $(argstr[1]): $(argstr[2]))"
            else
                "(std::max)({$(join(map(x->GenerateCppCode(x, info), args), ", "))})"
            end
        elseif name == "min"
            if length(args) == 2
                argstr = map(x->GenerateCppCode(x, info), args)
                "($(argstr[1]) < $(argstr[2])? $(argstr[1]): $(argstr[2]))"
            else
                "(std::min)({$(join(map(x->GenerateCppCode(x, info), args), ", "))})"
            end
        # unary math functions defined in jl2cpp.hpp
        elseif name in ["abs", "sqrt", "log", "log10", "exp", "sin", "cos", "tan", "asin", "acos", "atan", "sinh", "cosh", "tanh", "floor", "ceil", "round", "isnan", "isfinite", "isinf", "erf", "erfc"]
            "jl2cpp::$name($(GenerateCppCode(args[1], info)))"
        # cast
        elseif IsType(Symbol(name))
            "static_cast<$(Symbol(name)|>eval|>ConvertType)>($(GenerateCppCode(args[1], info)))"
        # just call the function if it is known user-defined function
        elseif name in info.userDefFunc
            "$name($(join(map(x->"jl2cpp::eval_if_possible(" * GenerateCppCode(x, info) * ")", args), ", ")))"
        else
            error("unknown function $name")
        end
    end

    """
    Convert Julia code into C++ code.
    """
    function GenerateCppCode(ex::Expr, info::GenCppCodeInfo)::String
        # operator or function call
        if ex.head == :call
            # assign operators which equivalent form and function to operators of C++
            if ex.args[1] in [:+, :-, :*, :/, :%, :|, :&, :<<, :<, :>, :(<=), :(>=), :(==), :(!=)]
                if length(ex.args) < 3
                    "($(ex.args[1])$(GenerateCppCode(ex.args[2], info)))"
                else
                    ret = foldl((x,y)->begin
                        y_ = GenerateCppCode(y, info)
                        "$x $(ex.args[1]) $y_"
                    end, GenerateCppCode(ex.args[2], info), ex.args[3:end])
                    "($ret)"
                end
                # convert operators which only Julia has
            elseif ex.args[1] in [:^, :$, :⊻, :.+, :.-, :.*, :./, :.^, :.<, :.<=, :.>, :.>=, :.==, :.!=, :>>, :>>>]
                arg1 = GenerateCppCode(ex.args[2], info)
                arg2 = GenerateCppCode(ex.args[3], info)
                if ex.args[1] == :^
                    "jl2cpp::pow($arg1, $arg2)" # TODO: implement matrix ^ scalar case (which is one of unique feature of julia)
                elseif ex.args[1] in [:&, :⊻] # old and new form of binary operator xor
                    "($arg1 ^ $arg2)"
                elseif ex.args[1] == :.+
                    "jl2cpp::add($arg1, $arg2)"
                elseif ex.args[1] == :.-
                    "jl2cpp::sub($arg1, $arg2)"
                elseif ex.args[1] == :.*
                    "jl2cpp::cwiseProduct($arg1, $arg2)"
                elseif ex.args[1] == :./
                    "jl2cpp::cwiseQuotient($arg1, $arg2)"
                elseif ex.args[1] == :.^
                    "jl2cpp::cwisePow($arg1, $arg2)"
                elseif ex.args[1] == :.<
                    "jl2cpp::cwiseLt($arg1, $arg2)"
                elseif ex.args[1] == :.<=
                    "jl2cpp::cwiseLte($arg1, $arg2)"
                elseif ex.args[1] == :.>
                    "jl2cpp::cwiseGt($arg1, $arg2)"
                elseif ex.args[1] == :.>=
                    "jl2cpp::cwiseGte($arg1, $arg2)"
                elseif ex.args[1] == :.==
                    "jl2cpp::cwiseEqual($arg1, $arg2)"
                elseif ex.args[1] == :.!=
                    "jl2cpp::cwiseNotEqual($arg1, $arg2)"
                elseif ex.args[1] == :>> # arith bitshift
                    "jl2cpp::arith_rshift($arg1, $arg2)"
                elseif ex.args[1] == :>>> # logical bitshift
                    "jl2cpp::logic_rshift($arg1, $arg2)"
                end
            else
                # TODO: make unknown and not included function binded with the parent code
                ConvertFunction("$(ex.args[1])", length(ex.args) > 1? ex.args[2:end] : [], info)
            end

        # assign
        # TODO: consider distinguish lhv and rhv

        elseif ex.head ==:(=)
            decl = ""
            ret = []
            rhs = GenerateCppCode(ex.args[end], info)
            i = -1
            # tuple in lhs
            if isa(ex.args[1], Expr) && ex.args[1].head == :tuple
                # multiple assignment
                if ex.args[2] isa Expr && ex.args[2].head == :tuple
                    # ad-hoc for special case: swap
                    if length(ex.args[1].args) == 2 && symdiff(ex.args[1].args, ex.args[2].args) |> isempty
                        return "std::swap(" * join(map(ex.args[1].args) do x GenerateCppCode(x, info); end, ", ") * ")"
                    end
                    for i = 1:(ex.args[1].args|>length)
                        if ex.args[1].args[i] isa Expr && ex.args[1].args[i].head == :local
                            push!(ret, "auto $(GenerateCppCode(ex.args[1].args[i].args[1], info)) = $(GenerateCppCode(ex.args[2].args[i], info))")
                        else
                            push!(ret, "$(GenerateCppCode(ex.args[1].args[i], info)) = $(GenerateCppCode(ex.args[2].args[i], info))")
                        end
                    end
                else
                    varname = []
                    for i=1:length(ex.args[1].args)
                        if ex.args[1].args[i] isa Expr && ex.args[1].args[i].head == :local
                            decl *= "typename jl2cpp::ElType<decltype(std::get<$(i - 1)>($rhs))>::type $(ex.args[1].args[i].args[1]); "
                            push!(varname, ex.args[1].args[i].args[1])
                        else
                            push!(varname, GenerateCppCode(ex.args[1].args[i], info))
                        end
                    end
                    push!(ret, "std::tie($(join(varname, ", "))) = $rhs")
                end
            else
                if ex.args[2] isa Expr && ex.args[2].head in [:local, :(=), :+=, :-=, :*=, :/=, :%=, :|=, :&=, :<<=, :>>=, :^=, :$=, :⊻=, :.+=, :.-=, :.*=, :./=, :.^=]
                    vn = if ex.args[2].head == :local ex.args[2].args[1].args[1] else ex.args[2].args[1] end
                    push!(ret, "$(GenerateCppCode(ex.args[2], info))")
                    push!(ret, "$(GenerateCppCode(ex.args[1], info)) = $vn")
                #if ex.args[1] isa Expr && ex.args[1].head == :local
                    # normal vardecl w/ init
                #    push!(ret, "auto $(GenerateCppCode(ex.args[1].args[1], info)) = $(GenerateCppCode(ex.args[2], info))")
                else
                    push!(ret, "$(GenerateCppCode(ex.args[1], info)) = $(GenerateCppCode(ex.args[2], info))")
                end
            end
            return decl * join(ret, "; ")# * (length(ret) <= 1? ";": "")

        # assign operators which equivalent form and function to operators of C++
        elseif ex.head in [:+=, :-=, :*=, :/=, :%=, :|=, :&=, :<<=, :>>=]
            "$(GenerateCppCode(ex.args[1], info)) $(ex.head) $(GenerateCppCode(ex.args[2], info))"
        elseif ex.head in [:^=, :$=, :⊻=, :.+=, :.-=, :.*=, :./=, :.^=]
            arg1 = GenerateCppCode(ex.args[1], info)
            arg2 = GenerateCppCode(ex.args[2], info)
            if ex.args[1] == :^
                "$arg1 = jl2cpp::pow($arg1, $arg2)"
            elseif ex.args[1] in [:&, :⊻] # old and new form of binary operator xor
                "$arg1 = $arg1 ^ $arg2"
            elseif ex.args[1] == :.+=
                "$arg1 = $arg1 + $arg2"
            elseif ex.args[1] == :.-=
                "$arg1 = $arg1 - $arg2"
            elseif ex.args[1] == :.*=
                "$arg1 = cwiseProduct($arg1, $arg2)"
            elseif ex.args[1] == :./=
                "$arg1 = cwiseQuotient($arg1, $arg2)"
            elseif ex.args[1] == :.^=
                "$arg1 = cwisePow($arg1, $arg2)"
            end
        # chained comparison operators
        # FixMe: operand with side effects causes wrong result
        elseif ex.head == :comparison
            ret = []
            for i = 2:2:length(ex.args)
                push!(ret, GenerateCppCode(ex.args[i-1]) * " " * "$(ex.args[i])" * " " * GenerateCppCode(ex.args[i+1]))
            end
            join(ret, " && ")
        # block
        elseif ex.head == :block
            ret = ""
            lineStr = ""
            info.indent += 1
            indentStr = repeat("\t", info.indent)
            for i = 1:length(ex.args)
                if isa(ex.args[i], Expr) && ex.args[i].head == :line
                    lineStr = GenerateCppCode(ex.args[i], info)
                    continue
                end
                l = GenerateCppCode(ex.args[i], info)
                ret *= l[end] == '}'? "$indentStr$lineStr\n$indentStr$l\n": "$indentStr$l; $lineStr\n"
                lineStr = ""
            end
            info.indent -= 1
            indentStr = repeat("\t", info.indent)
            "{\n$ret$indentStr}"
        elseif ex.head == :let
            body = GenerateCppCode(ex.args[1], info)[3:end]
            vars = length(ex.args) > 1? foldl((x,y)->"$x$(GenerateCppCode(y, info));\n", "", ex.args[2:end]): ""
            "{\n$vars$body"
        elseif ex.head == :(::)
            "$(ConvertType(ex.args[2]|>eval)) $(ex.args[1])"
        # function
        elseif ex.head == :function
            if ex.args[1].head != :(::)
                error("function must be typed")
            end
            name = ex.args[1].args[1].args[1]
            args = ex.args[1].args[1].args[2:end]
            retty = ex.args[1].args[2] |> eval |>ConvertType
            # retty = (retty.name.name == :Array)? "jl2cpp::matrix": ConvertType(retty)

            push!(info.userDefFunc, String(name))
            argstrs = String[]
            initstrs = String[]
            templatestrs = String[]
            i = 0
            for a = args
                if typeof(a) != Expr || a.head != :(::)
                    error("function arguments must be typed")
                end
                cname = a.args[1]
                ty = eval(a.args[2])
                templatename = "T$i"
                if ty.name.name == :Array
                    cty = ty.parameters[1] |> ConvertType
                    if ty.parameters[2] == 1
                        push!(argstrs, "$templatename& $(cname)")
                        #push!(initstrs, repeat("\t", info.indent + 1) * "Map<Matrix<$cty, Dynamic, 1> > $cname(($cty *)$(cname)_jl2cpp_arg->ptr, $(cname)_jl2cpp_arg->rsize);")
                        # TODO : uncomment
                        push!(initstrs, repeat("\t", info.indent + 1) * "$cty *const __restrict $(cname)_jl2cpp_ptr = ($cty *)$(cname).data(); auto const $(cname)_jl2cpp_rsize = $(cname).size();")
                        push!(info.argArrays, cname)
                    elseif ty.parameters[2] == 2
                        push!(argstrs, "$templatename& $(cname)")
                        #push!(initstrs, repeat("\t", info.indent + 1) * "Map<Matrix<$cty, Dynamic, Dynamic> > $cname(($cty *)$(cname)_jl2cpp_arg->ptr, $(cname)_jl2cpp_arg->rsize, $(cname)_jl2cpp_arg->csize);")
                        # TODO: uncomment
                        push!(initstrs, repeat("\t", info.indent + 1) * "$cty *const __restrict $(cname)_jl2cpp_ptr = ($cty *)$(cname).data(); auto const $(cname)_jl2cpp_rsize = $(cname).rows(); auto const $(cname)_jl2cpp_csize = $(cname).cols();")
                        push!(info.argArrays, cname)
                    else
                        error("3+ dimension array is currently not supported.")
                    end
                    push!(templatestrs, "typename $templatename")
                    i += 1
                else
                    cty = ty |> ConvertType
                    push!(argstrs, "$cty $cname")
                end
            end
            argstr = join(argstrs, ", ")
            initstr = join(initstrs, "\n")
            body = "{\n" * initstr * GenerateCppCode(ex.args[2], info)[2:end]
            (templatestrs|>isempty? "": "template<$(join(templatestrs, ", "))> ") * "EIGEN_STRONG_INLINE $retty $name($argstr)\n$body\n$(GenerateWrapperFunction(ex, info))"
        # for
        elseif ex.head == :for
            function makeFor(loopvar, iteratee)
                # foreach-style for
                if isa(iteratee, Symbol)
                    "for(auto&& $loopvar : $iteratee)"
                elseif isa(iteratee, Expr)
                    # number-counting for
                    if iteratee.head == :(:)
                        if length(iteratee.args) == 2
                            b = GenerateCppCode(iteratee.args[1], info)
                            e = GenerateCppCode(iteratee.args[2], info)
                            "for(auto $loopvar = $b; $loopvar <= $e; ++$loopvar)"
                        elseif length(iteratee.args) == 3
                            b = GenerateCppCode(iteratee.args[1], info)
                            s = GenerateCppCode(iteratee.args[2], info)
                            e = GenerateCppCode(iteratee.args[3], info)
                            if iteratee.args[2] isa Integer
                                if iteratee.args[2] > 0
                                    "for(auto $loopvar = $b; $loopvar <= $e; $loopvar += $s)"
                                else
                                    "for(auto $loopvar = $b; $loopvar >= $e; $loopvar += $s)"
                                end
                            else
                                "for(auto $loopvar = $b; ($s > 0)? ($loopvar <= $e): ($loopvar >= $e); $loopvar += $s)"
                            end
                        else
                            error("counting for with illegal range operator")
                        end
                    else
                        "for(auto&& $loopvar : $(GenerateCppCode(iteratee, info)))\n$(GenerateCppCode(loopbody, info))"
                    end
                end
            end
            if ex.args[1].head == :(=) # non-nested loop
                header = makeFor(ex.args[1].args[1], ex.args[1].args[2])
                body = GenerateCppCode(ex.args[2], info)
                "$header$body"
            elseif ex.args[1].head ==:block # nested loop
                header = foldl((x,y)-> x * repeat("\t", info.indent) * y, map(x->makeFor(x.args[1], x.args[2]), ex.args[1].args[1:end]))
                body = GenerateCppCode(ex.args[2], info)
                "$header$body"
            end

        elseif ex.head == :while
            cond = GenerateCppCode(ex.args[1], info)
            body = GenerateCppCode(ex.args[2], info)
            "while($cond)$body"
        # if and ternary operator
        elseif ex.head == :if
            # syntax if
            if isa(ex.args[2], Expr) && ex.args[2].head == :block
                cond = GenerateCppCode(ex.args[1], info)
                tclause = GenerateCppCode(ex.args[2], info)
                if length(ex.args) > 2
                    fclause = GenerateCppCode(ex.args[3], info)
                    "if($cond)$(tclause)else $fclause"
                else
                    "if($cond)$(tclause)"
                end
            else
                # ternary operator
                if length(ex.args) != 3
                    error("ternary conditional operator without false clause is not supported")
                end
                cond = GenerateCppCode(ex.args[1], info)
                tclause = GenerateCppCode(ex.args[2], info)
                fclause = GenerateCppCode(ex.args[3], info)

                "(($cond)? ($tclause): ($fclause))"
            end

        # array reference
        elseif ex.head == :ref
            isIdxScalar(idx) = typeof(idx) != Expr || idx.head != :(:) # scalar index value
            isIdxRange(idx) = !isIdxScalar(idx) && length(idx.args) == 3 # [a:b]
            isIdxRangeHead(idx) = isIdxRange(idx) && idx.args[1] == 1
            isIdxRangeTail(idx) = isIdxRange(idx) && idx.args[2] == :end
            isIdxFullRange(idx) = !isIdxScalar(idx) && (length(idx.args) == 1 || length(idx.args) == 2 && idx.args[1] == 1 && idx.args[2] == :end) # [:], [1:end]
            isIdxStride(idx) = !isIdxScalar(idx) && length(idx.args) == 3

            arrayName = GenerateCppCode(ex.args[1], info)

            # reference indeces
            # [idx] or [idxb:idxe] or [rb:re, cb:ce] or [r, c]
            # defined not as constant values but as functions because intended to do lazy evaluation
            rb(e) = GenerateCppCode(e.args[2].args[1], info)
            re(e) = GenerateCppCode(e.args[2].args[2], info)
            cb(e) = GenerateCppCode(e.args[3].args[1], info)
            ce(e) = GenerateCppCode(e.args[3].args[2], info)
            r(e) = GenerateCppCode(e.args[2], info)
            c(e) = GenerateCppCode(e.args[3], info)
            idx(e) = GenerateCppCode(e.args[2], info)
            idxb(e) = GenerateCppCode(e.args[2].args[1], info)
            idxe(e) = GenerateCppCode(e.args[2].args[2], info)

            if length(ex.args) == 1
                "$arrayName(0)"
            elseif length(ex.args) == 2 # one-dimension access
                if isIdxScalar(ex.args[2])                              # [i]
                    if Symbol(arrayName) in info.argArrays
                        "$(arrayName)_jl2cpp_ptr[$(idx(ex)) - 1]"
                    else
                        "$(arrayName).data()[$(idx(ex)) - 1]"
                    end
                elseif isIdxFullRange(ex.args[2])                       # [:]
                    "$arrayName"
                elseif isIdxRangeHead(ex.args[2])                       # [1:i]
                    "$arrayName.head($(idx(ex)))"
                elseif isIdxRangeTail(ex.args[2])                       # [i:end]
                    "$arrayName.tail($arrayName.size() - $(idxb(ex)) + 1)"
                elseif isIdxRange(ex.args[2])                           # [i:j]
                    "$arrayName.segment($(idxb(ex)) - 1, $(idxe(ex)) - $(idxb(idx)) + 1)"
                elseif isIdxStride(ex.args[2])
                    error("stride access is currently not supported.")
                else
                    error("unknown reference to array")
                end
            elseif length(ex.args) == 3 # two-dimension access
                if isIdxStride(ex.args[2]) || isIdxStride(ex.args[3])
                    error("stride access is currently not supported.")
                elseif isIdxFullRange(ex.args[2])
                    if isIdxFullRange(ex.args[3])                       # [:, :]
                        "$arrayName"
                    elseif isIdxRangeHead(ex.args[3])                   # [:, 1:b]
                        "$arrayName.leftCols($(ce(ex)))"
                    elseif isIdxRangeTail(ex.args[3])                   # [:, a:end]
                        "$arrayName.rightCols($arrayName.cols() - $(cb(ex)) + 1)"
                    elseif isIdxRange(ex.args[3])                       # [:, a:b]
                        "$arrayName.middleCols($(eb(ex)) - 1, $(ce(ex)) - $(cb(ex)) + 1)"
                    elseif isIdxScalar(ex.args[3])                      # [:, c]
                        "$arrayName.col($(c(ex)) - 1)"
                    else
                        error("this line must be unreachable.")
                    end
                elseif isIdxRangeHead(ex.args[2])
                    if isIdxFullRange(ex.args[3])                       # [1:b, :]
                        "$arrayName.topRows($(re(ex)))"
                    elseif isIdxRangeHead(ex.args[3])                   # [1:b, 1:d]
                        "$arrayName.topLeftCorner($(re(ex)), $(ce(ex)))"
                    elseif isIdxRangeTail(ex.args[3])                   # [1:b, c:end]
                        "$arrayName.topRightCorner($(rb(ex)), $arrayName.cols() - $(cb(ex)) + 1)"
                    elseif isIdxRange(ex.args[3])                       # [1:b, c:d]
                        "$arrayName.block(0, $(cb(ex)) - 1, $(re(ex)), $(ce(ex)) - $(cb(ex)) + 1)"
                    elseif isIdxScalar(ex.args[3])                      # [1:b, c]
                        "$arrayName.col($(c(ex)) - 1).head($(re(ex)))"
                    else
                        error("this line must be unreachable.")
                    end
                elseif isIdxRangeTail(ex.args[2])
                    if isIdxFullRange(ex.args[3])                       # [a:end, :]
                        "$arrayName.bottomRows($arrayName.rows() - $(rb(ex)) + 1)"
                    elseif isIdxRangeHead(ex.args[3])                   # [a:end, 1:d]
                        "$arrayName.bottomLeftCorner($arrayName.rows() - $(rb(ex)) + 1, $(ce(ex)))"
                    elseif isIdxRangeTail(ex.args[3])                   # [a:end, c:end]
                        "$arrayName.bottomRightCorner($arrayName.rows() - $(rb(ex)) + 1, $arrayName.cols() - $(cb(ex)) + 1)"
                    elseif isIdxRange(ex.args[3])                       # [a:end, c:d]
                        "$arrayName.block($(rb(ex)) - 1, $(ex.args[3].args[1]|> GenerateCppCode) - 1, $arrayName.rows() - $(rb(ex)) + 1, $(ce(ex)) - $(cb(ex)) + 1)"
                    elseif isIdxScalar(ex.args[3])                      # [a:end, c]
                        "$arrayName.col($(c(ex)) - 1).tail($arrayName.rows() - $(rb(ex)) + 1))"
                    else
                        error("this line must be unreachable.")
                    end
                elseif isIdxRange(ex.args[2])
                    if isIdxFullRange(ex.args[3])                       # [a:b, :]
                        "$arrayName.middleRows($(rb(ex)) - 1, $(re(ex)) - $(rb(ex)) + 1)"
                    elseif isIdxRangeHead(ex.args[3])                   # [a:b, 1:d]
                        "$arrayName.block($(rb(ex)) - 1, 0, $(re(ex)) - $(rb(ex)) + 1, $(ce(ex)))"
                    elseif isIdxRangeTail(ex.args[3])                   # [a:b, c:end]
                        "$arrayName.block($(rb(ex)) - 1, $(cb(ex)) - 1, $(re(ex)) - $(rb(ex)) + 1, $arrayName.cols() - $(cb(ex)) + 1)"
                    elseif isIdxRange(ex.args[3])                       # [a:b, c:d]
                        "$arrayName.block($(rb(ex)) - 1, $(cb(ex)) - 1, $(re(ex)) - $(rb(ex)) + 1, $(ce(ex)) - $(cb(ex)) + 1)"
                    elseif isIdxScalar(ex.args[3])                      # [a:b, c]
                        "$arrayName.col($(c(ex)) - 1).segment($(rb(ex)) - 1, $(re(ex)) - $(rb(ex)) + 1)"
                    else
                        error("this line must be unreachable.")
                    end
                elseif isIdxScalar(ex.args[2])
                    if isIdxFullRange(ex.args[3])                       # [a, :]
                        "$arrayName.row($(r(ex)) - 1)"
                    elseif isIdxRangeHead(ex.args[3])                   # [a, 1:d]
                        "$arrayName.row($(r(ex)) - 1).head($(cb(ex)))"
                    elseif isIdxRangeTail(ex.args[3])                   # [a, c:end]
                        "$arrayName.row($(r(ex)) - 1).tail($arrayName.cols() - $(cb(ex)) + 1)"
                    elseif isIdxRange(ex.args[3])                       # [a, c:d]
                        "$arrayName.row($(r(ex)) - 1).segment($(cb(ex)) - 1, $(ce(ex)) - $(cb(ex)) + 1)"
                    elseif isIdxScalar(ex.args[3])                      # [a, c]
                        if Symbol(arrayName) in info.argArrays
                            "$(arrayName)_jl2cpp_ptr[($(r(ex)) - 1) + ($(c(ex)) - 1) * $(arrayName)_jl2cpp_rsize]"
                        else
                            "$(arrayName).data()[($(r(ex)) - 1) + ($(c(ex)) - 1) * $(arrayName)_jl2cpp_rsize]"
                        end
                    end
                else
                    error("this line must be unreachable.")
                end
            else
                error("3+ dimeinson array is currently not supported.")
            end
            #procArray(arrayName, reverse(ex.args[2:end]))

        elseif ex.head == :local
            ret = []
            function getIval(ex)
                if ex.args[2] isa Expr && ex.args[2].head in [:local, :(=), :+=, :-=, :*=, :/=, :%=, :|=, :&=, :<<=, :>>=, :^=, :$=, :⊻=, :.+=, :.-=, :.*=, :./=, :.^=]
                    vn = if ex.args[2].head == :local ex.args[2].args[1].args[1] else ex.args[2].args[1] end
                    push!(ret, "$(GenerateCppCode(ex.args[2], info)) ")
                    "$vn"
                else
                    GenerateCppCode(ex.args[2], info)
                end
            end
            for i in 1:length(ex.args)
                if isa(ex.args[i], Expr)
                    if ex.args[i].head == :(=)
                        if isa(ex.args[i].args[1], Expr) && ex.args[i].args[1].head == :(::)
                            tname = ex.args[i].args[1].args[2] |> eval |> ConvertType
                            vname = ex.args[i].args[1].args[1]

                            ival = getIval(ex.args[i])
                            push!(ret, "$tname $vname = $ival")
                        else
                            vname = ex.args[i].args[1]
                            ival = getIval(ex.args[i])
                            push!(ret, "auto $vname = $ival")
                        end
                    elseif ex.args[i].head == :(::)
                        tname = ex.args[i].args[2] |> eval |> ConvertType
                        vname = ex.args[i].args[1]
                        push!(ret, "$tname $vname")
                    end
                else
                    # push!(ret, "auto $(ex.args[i])")
                    error("Local variable $(ex.args[i]) is declared without type annotation nor initialization.\n Please add at least one of them.")
                end
            end
            join(ret, "; ")
        elseif ex.head == :tuple
            # a,b,c=... form not in :block
            if isa(ex.args[end], Expr) && ex.args[end].head == :(=)
                decl = ""
                tie = "std::tie("
                rhs = GenerateCppCode(ex.args[end].args[end], info)
                for i = 1:(ex.args|>length) - 1
                    if isa(ex.args[i], Expr) && ex.args[i].head == :local
                        decl *= "decltype(std::get<$i>($rhs) $(GenerateCppCode(ex.args[i].args[1], info)); "
                        tie *= GenerateCppCode(ex.args[i].args[1], info)
                    else
                        tie *= GenerateCppCode(ex.args[i], info)
                    end
                    tie *= ", "
                end
                if isa(ex.args[end].args[end], Expr) && ex.args[end].args[1].head == :local
                    decl *= "decltype(std::get<$(ex.args|>length)>($rhs)) $(GenerateCppCode(ex.args[end].args[1].args[1], info)); "
                    tie *= GenerateCppCode(ex.args[end].args[1].args[1], info)
                else
                    tie *= GenerateCppCode(ex.args[end].args[1], info)
                end
                "$decl$tie) = $rhs"
            else
                # rhs tuple
                # "std::forward_as_tuple(" * join(map(x->GenerateCppCode(x, info), ex.args), ", ") * ")"
                "std::make_tuple(" * join(map(x->GenerateCppCode(x, info), ex.args), ", ") * ")"
            end
        # vector constructor
        # TODO: correct :vect, :hcat, :vcat
        elseif ex.head == :vcat
            if isa(ex.args[1], Expr) && ex.args[1].head == :row #matrix
                el = []
                for c in ex.args, r in c.args
                    push!(el, GenerateCppCode(r, info))
                end
                "(Matrix<jl2cpp::ElType<decltype(jl2cpp::get_one_elem($(el[1])))>::type, Dynamic, Dynamic>(jl2cpp::get_rows($(el[1])) * $(length(ex.args)),  jl2cpp::get_cols($(el[1])) * $(length(ex.args[1].args))) << $(join(el, ", "))).finished()"
            else # vertical vector
                el = []
                for c in ex.args
                    if c isa Expr && c.head == :(:) # TODO: implement range expression in vcat
                        error("range expression is currently not supported")
                    end
                    push!(el, GenerateCppCode(c, info))
                end
                "(Matrix<jl2cpp::ElType<decltype(jl2cpp::get_one_elem($(el[1])))>::type, Dynamic, Dynamic>(jl2cpp::get_rows($(el[1])) * $(length(ex.args)), jl2cpp::get_cols($(el[1]))) << $(join(el, ", "))).finished()"
            end
        elseif ex.head == :hcat # horizontal vector
            if ex.args|>isempty
                # empty array
                error("not implemented")
            else
                el = []
                for c in ex.args
                    push!(el, GenerateCppCode(c, info))
                end
                "(Matrix<jl2cpp::ElType<decltype(jl2cpp::get_one_elem($(el[1])))>::type, Dynamic, Dynamic>(jl2cpp::get_rows($(el[1])), jl2cpp::get_cols($(el[1])) * $(length(ex.args))) << $(join(el, ", "))).finished()"
            end
        elseif ex.head == :vect #1d array
            if ex.args|>isempty
                # empty array
                error("not implemented")
            else
                el = []
                for c in ex.args
                    push!(el, GenerateCppCode(c, info))
                end
                "(Matrix<jl2cpp::ElType<decltype($(el[1]))>::type, Dynamic, 1>($(length(ex.args))) << $(join(el, ", "))).finished()"
            end
        elseif ex.head == Symbol(".'")
            "$(GenerateCppCode(ex.args[1], info)).transpose()"
        elseif ex.head == Symbol("'")
            "$(GenerateCppCode(ex.args[1], info)).adjoint()"
        # macro
        elseif ex.head == :macrocall
            if ex.args[1] == Symbol("@fastmath")
                GenerateCppCode(ex.args[2], info)
            elseif ex.args[1] == Symbol("@printf")
                "printf($(join(map(x->GenerateCppCode(x, info), ex.args[2:end]), ", ")))"
            elseif ex.args[1] == Symbol("@test")
                "assert($(join(map(x->GenerateCppCode(x, info), ex.args[2:end]), " && ")))"
            elseif ex.args[1] == Symbol("@pragma")
                "#pragma $(ex.args[2])\n" * repeat("\t", info.indent) * GenerateCppCode(ex.args[3], info)
            elseif ex.args[1] == Symbol("@parallel")
                if !(ex.args[2] isa Expr && ex.args[2].head == :for)
                    error("not implemented")
                end
            elseif ex.args[1] == (Symbol("@threads"))
                "#pragma omp parallel for\n" * repeat("\t", info.indent) * GenerateCppCode(ex.args[2], info)
            elseif ex.args[1] isa Expr && ex.args[1].head == :(.) # module
                if ex.args[1].args[1] == :Threads && ex.args[1].args[2] == Expr(:quote, Symbol("@threads"))
                    "#pragma omp parallel for\n" * repeat("\t", info.indent) * GenerateCppCode(ex.args[2], info)
                else
                    "/*$(ex.args[1])*/ " * GenerateCppCode(ex.args[2], info)
                end
            else "/*$(ex.args[1])*/ " * GenerateCppCode(ex.args[2], info)
            end
        elseif ex.head == :using
        # line number for debug
        elseif ex.head == :line
            "/* line $(ex.args[1]) */"
        elseif ex.head == :return
            "return jl2cpp::eval_if_possible($(GenerateCppCode(ex.args[1], info)))"
        elseif ex.head == :continue
            "continue"
        elseif ex.head == :break
            "break"
        else
            "/** unknown expr: $(ex.head) **/"
        end
    end
    function GenerateCppCode(e::AbstractString, info::GenCppCodeInfo)::String
        return "\"" * e * "\""
    end
    function GenerateCppCode(e::Any, info::GenCppCodeInfo)::String
        return string(e)
    end

    # insert `return` to AST
    function InsertReturnSyntax(ex :: Expr)
        if ex.head == :function
            ex.args[2] = InsertReturnSyntax(ex.args[2])
            ex
        elseif ex.head == :block
            ex.args[end] = InsertReturnSyntax(ex.args[end])
            ex
        elseif ex.head == :let
            InsertReturnSyntax(ex.args[1])
            ex
        elseif ex.head == :if
            # syntax if
            if isa(ex.args[2], Expr)
                ex.args[2] = InsertReturnSyntax(ex.args[2])
                ex.args[3] = InsertReturnSyntax(ex.args[3])
            else
                # ternary operator
                ex = Expr(:return, ex)
            end
            ex
        elseif ex.head in [:while, :for, :return]
            ex
        else
            # if the cluause that should be returned is a assignment, discard lvalue and modify it to return syntax.
            # this may cause compatibility issue in (very rare) case that the lvalue is a global variable or causes side effect by the assignment.
            if isa(ex, Expr) && ex.head == :(=)
                Expr(:return, ex.args[2])
            else
                Expr(:return, ex)
            end
        end
    end

    function InsertReturnSyntax(ex)
        Expr(:return, ex)
    end

    function IsVarExist(env :: Vector{Dict{Symbol, DataType}}, var :: Symbol) :: Bool
        if env|>isempty
            return false
        end
        for i=length(env):-1:1
            if haskey(env[i], var)
                return true
            end
        end
        false
    end

    """
    Insert explicit `local` node to the first assignment of local variables
    in order to make variable decleration easy to be converted.
    """
    function InsertVarDecl(ex::Expr, env :: Vector{Dict{Symbol, DataType}})
        function NewScope!(env :: Vector{Dict{Symbol, DataType}})
            push!(env, Dict{Symbol, DataType}())
        end
        function DeleteScope!(env :: Vector{Dict{Symbol, DataType}})
            pop!(env)
        end
        function RegisterVar!(env :: Vector{Dict{Symbol, DataType}}, var::Any, typ::DataType)
            if isempty(env)
                NewScope!(env)
            end
            if typeof(var) != Symbol
                println(var)
                join(stacktrace(),"\n") |> println
                error("not symbol")
            end
            env[end][var] = typ
        end

        if ex.head == :function
            if isa(ex.args[1], Expr) && ex.args[1].head == :(::)
                NewScope!(env)
                for x in ex.args[1].args[1].args
                    if isa(x, Expr) && x.head == :(::)
                        RegisterVar!(env, x.args[1], x.args[2]|>eval)
                    else
                        RegisterVar!(env, x, Any)
                    end
                end
                InsertVarDecl(ex.args[2], env)
                DeleteScope!(env)
            else
                NewScope!(env)
                if length(ex.args[1].args) >= 2
                    # println(ex.args[1].args)
                    foldl((x,y)->RegisterVar!(env, y, Any), ex.args[1].args[1:end])
                end
                InsertVarDecl(ex.args[2], env)
                DeleteScope!(env)
            end

        elseif ex.head == :for
            NewScope!(env)
            if ex.args[1].head == :block
                foldl((x,y)->RegisterVar!(env, y.args[1], Any), ex.args[1].args[1:end])
            else
                RegisterVar!(env, ex.args[1].args[1], Any)
            end
            InsertVarDecl(ex.args[2], env)
            DeleteScope!(env)
        elseif ex.head == :while
            NewScope!(env)
            if ex.args[1].head == :block
                foldl((x,y)->RegisterVar!(env, y.args[1], Any), ex.args[1].args[1:end])
            else
                RegisterVar!(env, ex.args[1].args[1], Any)
            end
            InsertVarDecl(ex.args[2], env)
            DeleteScope!(env)
        elseif ex.head == :block
            for i in 1:length(ex.args)
                if ex.args[i].head == :line
                    continue
                else
                    ex.args[i] = InsertVarDecl(ex.args[i], env)
                end
            end
        elseif ex.head == :(=)
            ex.args[2] = InsertVarDecl(ex.args[2], env)
            if isa(ex.args[1], Expr) && ex.args[1].head == :(::) && isa(ex.args[1].args[1], Symbol)
                if !IsVarExist(env, ex.args[1].args[1])
                    RegisterVar!(env, ex.args[1].args[1], ex.args[1].args[2] |> eval)
                    ex = Expr(:local, ex)
                end
            elseif isa(ex.args[1], Expr) && ex.args[1].head == :tuple
                for i in 1:(ex.args[1].args|>length)
                    if !IsVarExist(env, ex.args[1].args[i])
                        RegisterVar!(env, ex.args[1].args[i], Any)
                        ex.args[1].args[i] = Expr(:local, ex.args[1].args[i])
                    end
                end
            else
                if isa(ex.args[1], Symbol) && !IsVarExist(env, ex.args[1])
                    RegisterVar!(env, ex.args[1], Any) # TODO : type inference
                    ex = Expr(:local, ex)
                end
            end
        elseif ex.head == :tuple && isa(ex.args[end], Expr) && ex.args[end].head == :(=)
            for i = 1 : ex.args |> length - 1
                if !IsVarExist(env, ex.args[i])
                    RegisterVar!(env, ex.args[i], Any)
                    ex = Expr(:local, ex.args[i])
                end
            end
            if !IsVarExist(env, ex.args[end].args[1])
                RegisterVar!(env, ex.args[end].args[1], Any)
                ex.args[end].args[1] = Expr(:local, ex.args[end].args[1])
            end
        elseif ex.head == :macrocall
            for i = 2:length(ex.args)
                InsertVarDecl(ex.args[i], env)
            end
        elseif ex.head == :local
            for i = 1:length(ex.args)
                if ex.args[i] isa Symbol
                    RegisterVar!(env, ex.args[i], Any)
                    continue
                elseif ex.args[i] isa Expr
                    e = ex.args[i]
                    if e.head == :(::)
                        RegisterVar!(env, e.args[1], e.args[2]|>eval)
                        continue
                    elseif e.head == :(=)
                        InsertVarDecl(e.args[2], env)
                        lhs = e.args[1]
                        if lhs isa Expr && lhs.head == :(::)
                            RegisterVar!(env, lhs.args[1], lhs.args[2]|>eval)
                            continue
                        elseif lhs isa Symbol
                            RegisterVar!(env, lhs, Any)
                            continue
                        else
                            error("this cannot be reached")
                        end
                    else
                        error("wrong syntax")
                    end
                else
                    error("wrong syntax")
                end
            end
        end
        ex
    end
    InsertVarDecl(ex :: Expr) = InsertVarDecl(ex, Vector{Dict{Symbol, DataType}}())
    InsertVarDecl(ex :: Any, env :: Vector{Dict{Symbol, DataType}}) = ex

    # make form of `f(x)=x*2` function to `function f(x); x*2; end`
    function makeShortSyntaxFunctionLong(ex :: Expr)
        if ex.head == :(=) && ex.args[1] isa Expr &&
            (ex.args[1].head == :call ||
            (ex.args[1].head == :(::) && ex.args[1].args[1] isa Expr && ex.args[1].args[1].head == :call)
            )
                ex.head = :function
        end
        map!(makeShortSyntaxFunctionLong, ex.args, ex.args)
        ex
    end
    makeShortSyntaxFunctionLong(x) = x

    # replace "!" in function name to "_"
    function sanitizeFuncName(ex::Expr)
        if ex.head == :call && ex.args[1] isa Symbol
            ex.args[1] = Symbol(replace(string(ex.args[1]), "!", "_"))
            if ex.args|>length > 2
                map!(sanitizeFuncName, ex.args[2:end], ex.args[2:end])
            end
        else
            map!(sanitizeFuncName, ex.args, ex.args)
        end
        ex
    end
    sanitizeFuncName(x) = x

    function WriteCppFile(src::String, filename::String)
        open(filename * ".cpp", "a") do os
            println(os, src)
        end
    end
    function CompileCpp(filename :: String)
        if haskey(ENV, "NO_COMPILE") && ENV["NO_COMPILE"] == "1"
            info("no compiling")
            return
        end
        ext = Libdl.dlext
        run(`$CC $C_OPT $filename.cpp -o $filename.$ext`)
    end

    #TODO: function translation
    #TODO: type inference
    function ConvertFunctionDefinition(ex :: Expr, filename::String, info::GenCppCodeInfo)
        #ex|>Meta.show_sexpr
        if ex.head != :function
            error("this cannot be reached")
        end
        if ex.args[1].head != :(::)
            error("functions must be typed")
        end

        origname = ex.args[1].args[1].args[1]
        funname = replace(string(ex.args[1].args[1].args[1]), "!", "_")
        funname = "$(funname)_jl2cpp_wrapper"
        retty = eval(ex.args[1].args[2])

        if retty.name.name != :Void
            ex = ex |> InsertReturnSyntax
        end

        code = GenerateCppCode(ex |> sanitizeFuncName |> InsertVarDecl, info)
        WriteCppFile(code, filename)

        callexpr = Expr(:call, :ccall,
            Expr(:tuple, Expr(:quote, Symbol(funname)), joinpath(pwd(), "$filename." * Libdl.dlext)),
            retty.name.name == :Array? Expr(:curly, :(Jl2Cpp.CMatrix), retty.parameters[1]): retty.name.name == :Tuple? :(Jl2Cpp.CTuple): retty,
            Expr(:tuple)
        )
        freeexpr = Expr(:call, :ccall,
        Expr(:tuple, Expr(:quote, Symbol(funname)), joinpath(pwd(), "$filename." * Libdl.dlext)),
        Void,
        Expr(:tuple)
    )
        # repair original function name because might be modified
        ex.args[1].args[1].args[1] = origname
        precall_code = Expr[]
        push!(precall_code, :(gc_enable(false)))
        for arg in ex.args[1].args[1].args[2:end]
            if arg.head != :(::)
                error("all arguments must be typed")
            end

            ty = eval(arg.args[2])
            if ty.name.name == :Array
                conv_code = Expr(:(=))
                if ty.parameters[2] == 1
                    conv_code = :(tmp = Ref(Jl2Cpp.CMatrix{$(ty.parameters[1])}(pointer($(arg.args[1])), length($(arg.args[1])), 1)))
                    push!(callexpr.args[4].args, :(Ref{Jl2Cpp.CMatrix{$(ty.parameters[1])}}))
                elseif ty.parameters[2] == 2
                    conv_code = :(tmp = Ref(Jl2Cpp.CMatrix{$(ty.parameters[1])}(pointer($(arg.args[1])), size($(arg.args[1]), 1), size($(arg.args[1]), 2))))
                    push!(callexpr.args[4].args, :(Ref{Jl2Cpp.CMatrix{$(ty.parameters[1])}}))
                else
                    error("3+ dimension array is currently not supported")
                end
                conv_code.args[1] = Symbol("$(arg.args[1])_jl2cpp_arg")
                push!(precall_code, conv_code)
                push!(callexpr.args, Symbol("$(arg.args[1])_jl2cpp_arg"))

                # push!(precall_code, :(println($(Symbol("$(arg.args[1]i)_jl2cpp_arg")))))
            else
                push!(callexpr.args[4].args, arg.args[2])
                push!(callexpr.args, arg.args[1])
            end
        end

        postcall_code = []
        if retty.name.name == :Array
            push!(postcall_code, :(jl2cpp_ret = Jl2Cpp.ConvertToArray(jl2cpp_ret)))
        elseif retty.name.name == :Tuple
            push!(postcall_code, :(jl2cpp_ret = Jl2Cpp.ConvertToTuple($retty, jl2cpp_ret)))
        end
        #push!(postcall_code, :(gc_enable(true)))

        ex.args[2] = Expr(:block)
        append!(ex.args[2].args, precall_code)
        if retty.name.name == :Void
            push!(ex.args[2].args, callexpr)
            append!(ex.args[2].args, postcall_code)
        else
            push!(ex.args[2].args, Expr(:(=), Symbol("jl2cpp_ret"), callexpr))
            append!(ex.args[2].args, postcall_code)
            push!(ex.args[2].args, Expr(:return, Symbol("jl2cpp_ret")))
        end
        ex
    end

    function CppMacroImpl(ex::Expr, filename::String, info::GenCppCodeInfo)
        if ex.head == :block
            for i in 1:length(ex.args)
                CppMacroImpl(ex.args[i], filename, info)
            end
        elseif ex.head == :function
            ex = ConvertFunctionDefinition(ex, filename, info)
        elseif ex.head == :macrocall
            if ex.args[1] == Symbol("@inline") # ignore @inline because inline prefix is always added in C++ code
                CppMacroImpl(ex.args[2], filename, info)
            end
        elseif ex.head == :(=)
            error("not implemented")
        # elseif ex.head == :line
        else
            return ex
        end
        return ex
    end

    macro cpp(ex::Expr)
        putInfo("CC: $CC")
        putInfo("C_OPT: $C_OPT")

        filename = "jl2cpp_"*(PROGRAM_FILE|>basename|>splitext)[1]

        if isfile(filename * ".cpp")
            rm(filename * ".cpp")
        end
        if haskey(ENV, "JL2CPP_DEBUG") || "-D" in ARGS
            println("Debug mode")
        else
            push!(CPP_PREDEFS, "#define EIGEN_NO_DEBUG 1")
        end
        # cppfree = DLL_EXPORT * "void free_ptr(void* p){free(p);}\n"
        #WriteCppFile(join(CPP_PREDEFS, "\n") * "\n" * join(CPP_INCLUDES, "\n") * "\n" * join(CPP_USINGS, "\n") * "\n" * cppfree, filename)
        WriteCppFile(join(CPP_PREDEFS, "\n") * "\n" * join(CPP_INCLUDES, "\n") * "\n" * join(CPP_USINGS, "\n") , filename)

        ex = ex|>makeShortSyntaxFunctionLong
        for f in RECIPES
            f(ex)
        end
        info = GenCppCodeInfo(0, Set([]), Set([]))
        ret = CppMacroImpl(ex, filename, info)
        CompileCpp(filename)
        esc(ret)
    end
    export @cpp

    """
    Make definitions of both original version of passed function and compiled version of that.
    """
    macro test_func(ex :: Expr)
        ret = Expr(:block, copy(ex))
        name=ex.args[1].args[1].args[1]
        ex.args[1].args[1].args[1] = Symbol(name, "_cpp")
        # push!(ret.args, :(println("Start compiling ", $name)))
        # push!(ret.args, :(flush(STDOUT)))
        # push!(ret.args, :(tic()))
        push!(ret.args, :(@Jl2Cpp.cpp $ex))
        # push!(ret.args, :(println("Done. time: $(toq())s")))
        # push!(ret.args, :(flush(STDOUT)))
        esc(ret)
    end
    export test_func

    macro test_run(times :: Int64, fun :: Expr, cpp_only=false)
        prepare_data = Expr(:block)
        orig_expr = Expr(:block)
        cpp_expr = Expr(:block)
        call_orig = Expr(:(=), :orig_ret, Expr(:call, fun.args[1]))
        call_cpp = Expr(:(=), :cpp_ret, Expr(:call, Symbol(fun.args[1], "_cpp")))
        valid_args = Expr(:block)
        i = 0
        for a in fun.args[2:end]
            push!(prepare_data.args, Expr(:(=), Symbol("arg", i), a))
            push!(orig_expr.args, Expr(:(=), Symbol("orig_arg", i), Expr(:call, :copy,Symbol("arg", i))))
            push!(cpp_expr.args, Expr(:(=), Symbol("cpp_arg", i), Expr(:call, :copy,Symbol("arg", i))))
            push!(call_orig.args[2].args, Symbol("orig_arg", i))
            push!(call_cpp.args[2].args, Symbol("cpp_arg", i))
            push!(valid_args.args,
                Expr(:if, Expr(:call, :!, Expr(:call, :isapprox, Symbol("cpp_arg", i), Symbol("orig_arg", i))),
                    Expr(:call, :warn, "argument #", i+1, " is not expected value.")
                )
            )
            i += 1
        end
        push!(orig_expr.args, call_orig)
        push!(cpp_expr.args, call_cpp)
        quote
            local orig_time = 0.0
            local cpp_time = 0.0
            for i=1:$times
                info("Start iteration #$i")
                flush(STDOUT)
                $prepare_data
                if !$cpp_only
                    info("Runninng original")
                    flush(STDOUT)
                    tic()
                    $orig_expr
                    orig_time += toq()
                end
                info("Runninng cpp ver")
                flush(STDOUT)
                tic()
                $cpp_expr
                cpp_time += toq()
                if !$cpp_only
                    $valid_args
                    if cpp_ret != orig_ret && !(cpp_ret ≈ orig_ret)
                        warn("return value is wrong.")
                    end
                end
            end
            if !$cpp_only
                orig_time /= $times
                cpp_time /= $times
                println("orig: $(orig_time)s, cpp: $(cpp_time)s")
                println("cpp/orig ratio = ", cpp_time/orig_time)
                (orig_time, cpp_time)
            else
                cpp_time /= $times
                println("cpp: $(cpp_time)s")
                (0.0, cpp_time)
            end
        end|>esc
    end
    export test_run
end

#=
TODO:

=#


text:
# 测试稳态与死亡率关系的脚本

t = env.T

@time EasySearch.SteadyState!( t, Guess,
   Dt, Dst, Pt, Ps, Pc, env,
   atol = 1E-8,  # tolerance of Gauss-Seidel iteration
   MaxIter = 50,  # maximum loops
   PrintMode = "final",  # mode of printing
   MagicNum = 2.0,  # magic number, the lower bound of K/L (capital per labor)
   StepLen = 0.5  # relative step length to update guesses, in range (0,1]
)
PyPlot.subplot(2,2,1)
   PyPlot.plot(Ps[:F][t,:])
PyPlot.subplot(2,2,2)
   PyPlot.plot(Dst[:𝒜][t,:])
PyPlot.subplot(2,2,3)
   PyPlot.plot(Ps[:N][t,:])
PyPlot.tight_layout()

# # 3. plotting & output
# EasyPlot.Plot_SteadyState( 1, Dt, Dst, Pt, Ps, Pc, env,
#    outpdf = string("./output/", "InitSS_", EasyIO.LogTag(), ".pdf" )
#    )

PyPlot.plot(diff(Ps[:N][env.T,:]) )

#


text:
using Test
using ProbNumDiffEq
using OrdinaryDiffEq
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

@testset "Test the different diffusion models" begin
    prob = prob_ode_fitzhughnagumo
    true_sol = solve(prob, Tsit5(), abstol=1e-12, reltol=1e-12)

    @testset "Time-Varying Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=DynamicDiffusion(), smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Fixed Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=FixedDiffusion(), smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Fixed Diffusion - uncalibrated and with custom initial value" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=FixedDiffusion(1e3, false), smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Varying Diagonal Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=DynamicMVDiffusion(), smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Fixed Diagonal Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=FixedMVDiffusion(), smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Fixed Diagonal Diffusion - uncalibrated and with custom values" begin
        d = length(prob.u0)
        initial_diffusion = 1 .+ rand(d)
        sol = solve(
            prob,
            EK0(diffusionmodel=FixedMVDiffusion(initial_diffusion, false), smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end
end


text:
mutable struct Tree{T<:Integer} <: AbstractAgent
  id::T
  pos::Tuple{T, T}
  status::Bool  # true is green and false is burning
end

# we can put the model initiation in a function
function model_initiation(;f, d, p, griddims, seed)
  Random.seed!(seed)

  space = Space(griddims, moore = true)

  properties = Dict(:f => f, :d => d, :p => p)
  forest = ABM(Tree, space; properties=properties, scheduler=random_activation)

  # create and add trees to each node with probability d, which determines the density of the forest
  for node in 1:nv(forest)
    pp = rand()
    if pp <= forest.properties[:d]
      tree = Tree(node, (1,1), true)
      add_agent!(tree, node, forest)
    end
  end
  return forest
end

function forest_step!(forest)
  for node in nodes(forest; by = :random)
    if length(forest.space.agent_positions[node]) == 0  # the cell is empty, maybe a tree grows here?
      p = rand()
      if p <= forest.properties[:p]
        bigest_id =  maximum(keys(forest.agents))
        treeid = bigest_id +1
        tree = Tree(treeid, (1,1), true)
        add_agent!(tree, node, forest)
      end
    else
      treeid = forest.space.agent_positions[node][1]  # id of the tree on this cell
      tree = id2agent(treeid, forest)  # the tree on this cell
      if tree.status == false  # if it is has been burning, remove it.
        kill_agent!(tree, forest)
      else
        f = rand()
        if f <= forest.properties[:f]  # the tree ignites on fire
          tree.status = false
        else  # if any neighbor is on fire, set this tree on fire too
          neighbor_cells = node_neighbors(tree, forest)
          for cell in neighbor_cells
            treeid = get_node_contents(cell, forest)
            if length(treeid) != 0  # the cell is not empty
              treen = id2agent(treeid[1], forest)
              if treen.status == false
                tree.status = false
                break
              end
            end
          end
        end
      end
    end
  end
end