text:
export TimeDelay, ps2μm!, μm2ps!

ps2μm(t::Real) = round(149.896225 * t)
μm2ps(d::Real) = d / 149.896225

for f in (:ps2μm, :μm2ps)
    @eval begin
        function $(Symbol(f, !))(arr::VecI)
            @simd for i in eachindex(arr)
                @inbounds arr[i] = $f(arr[i])
            end
            return nothing
        end
    end
end

struct TimeDelay <: TimeAxis
    dat::Vector{Float64} # unit: μm

    function TimeDelay(src::VecI, unit::String="ps")
        dat = similar(src)
        if unit ≡ "ps"
            @simd for i in eachindex(dat)
                @inbounds dat[i] = ps2μm(src[i])
            end
        elseif unit ≡ "μm"
            @simd for i in eachindex(dat)
                @inbounds dat[i] = src[i]
            end
        else
            error("TimeDelay(..., unit = $unit) is not allowed.")
        end
        return new(dat)
    end

    function TimeDelay(s::String, unit::String="μm"; Itdy::Int=1)
        raw = readdlm(s)
        dat = Vector{Float64}(undef, size(raw, 1))
        if unit ≡ "μm"
            @simd for i in eachindex(dat)
                @inbounds dat[i] = raw[i, Itdy]
            end
        elseif unit ≡ "ps"
            @simd for i in eachindex(dat)
                @inbounds dat[i] = ps2μm(raw[i, Itdy])
            end
        else
            error("TimeDelay(..., unit = $unit) is not allowed.")
        end
        return new(dat)
    end
end


text:
using Plots,StatsPlots,VegaLite
import Plots: cgrad

include("Color.jl")

iscircfactor(f::Symbol) = f in [:Ori,:Dir,:Ori_Final,:Dir_Final,:HueAngle,:Angle]
factorunit(fs::Vector{Symbol};timeunit=SecondPerUnit)=join(factorunit.(fs,timeunit=timeunit),", ")
function factorunit(f::Symbol;timeunit=SecondPerUnit)
    fu=String(f)
    if any(occursin.(["Ori","Angle"],fu))
        fu="$fu (deg)"
    elseif fu=="dir"
        fu="Direction (deg)"
    elseif fu=="Diameter"
        fu="$fu (deg)"
    elseif fu=="SpatialFreq"
        fu="$fu (cycle/deg)"
    elseif fu=="TemporalFreq"
        fu = "$fu (cycle/sec)"
    elseif fu=="Time"
        if timeunit ==1
            fu="$fu (s)"
        elseif timeunit == 0.001
            fu="$fu (ms)"
        end
    elseif fu=="Response"
        fu="$fu (spike/s)"
    elseif fu=="ResponseF"
        fu="Response (% \\DeltaF / F)"
    end
    return fu
end
factorunitstring(f,u) = factorunitstring(string(f),string(u))
function factorunitstring(f::String,u::String)
    isempty(f) && return f
    isempty(u) && return f
    "$f ($u)"
end

"scatter plot of spike trains"
function plotspiketrain(x,y;group::Vector=[],timeline=[0],color=huecolors(length(unique(group))),title="",size=(800,600))
    nt = isempty(x) ? 0 : maximum(y)
    ms =0.45*size[2]/(nt+5)
    p = plot(;size,leg=false,title,grid=false)
    if isempty(group)
        scatter!(p,x,y;label="SpikeTrain",markershape=:vline,markersize=ms,markerstrokewidth = 0,color=RGBA(0,0,0,0.7))
    else
        scatter!(p,x,y;group,markershape=:vline,markersize=ms,markerstrokewidth = 0,color=permutedims(color))
    end
    vline!(p,timeline;line=(:grey),label="TimeLine",xaxis=(factorunit(:Time)),yaxis=("Trial"))
end
function plotspiketrain(sts::Vector;uids::Vector=[],sortvalues=[],timeline=[0],color=huecolors(0),title="",size=(800,600))
    if isempty(uids)
        g=uids;uc=color
    else
        fuids = flatspiketrains(uids,sortvalues)[1]
        g=map(i->"U$i",fuids);uc=huecolors(length(unique(fuids)))
    end
    plotspiketrain(flatspiketrains(sts,sortvalues)[1:2]...;group=g,timeline,color=uc,title,size)
end

# function plotspiketrain1(x::Vector,y::Vector,c::Vector=[];xmin=minimum(x)-10,xmax=maximum(x)+10,xgroup::Vector=[],
#     ymin=minimum(y)-1,ymax=maximum(y)+1,timemark=[0],theme=Theme(),
#     colorkey="",colorfun=Scale.lab_gradient(colorant"white",colorant"red"),colorminv=[],colormaxv=[])
#     xl="Time (ms)";yl="Trial"
#     if isempty(c)
#         if isempty(xgroup)
#             plot(x=x,y=y,xintercept=timemark,theme,Geom.point,Geom.vline(color="gray",size=1pt),
#             Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),Guide.xlabel(xl),Guide.ylabel(yl))
#         else
#             plot(x=x,y=y,xgroup=xgroup,xintercept=fill(timemark[1],length(x)),theme,
#             Geom.subplot_grid(Geom.point,Geom.vline(color="gray",size=1pt),free_x_axis=true),
#             Guide.xlabel(xl),Guide.ylabel(yl))
#         end
#     else
#         yl="$yl Sorted"
#         if isempty(colorminv);colorminv=minimum(c);end
#         if isempty(colormaxv);colormaxv=maximum(c);end
#         if isempty(xgroup)
#             plot(x=x,y=y,color=c,xintercept=timemark,theme,Geom.point,Geom.vline(color="gray",size=1pt),
#             Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),Guide.xlabel(xl),Guide.ylabel(yl),Guide.colorkey(colorkey),
#             Scale.ContinuousColorScale(colorfun,minvalue=colorminv,maxvalue=colormaxv))
#         else
#             plot(x=x,y=y,color=c,xgroup=xgroup,xintercept=fill(timemark[1],length(x)),theme,
#             Geom.subplot_grid(Geom.point,Geom.vline(color="gray",size=1pt),free_x_axis=true),
#             Guide.xlabel(xl),Guide.ylabel(yl),Guide.colorkey(colorkey),
#             Scale.ContinuousColorScale(colorfun,minvalue=colorminv,maxvalue=colormaxv))
#         end
#     end
# end
# plotspiketrain1(rvs::RVVector;sortvar=[],xgroup::Vector=[],timemark=[0],theme=Theme(),colorkey="",colorfun=Scale.lab_gradient(colorant"white",colorant"red"),colorminv=[],colormaxv=[]) = plotspiketrain1(flatrvs(rvs,sortvar)...,xgroup=xgroup,timemark=timemark,theme=theme,colorkey=colorkey,colorfun=colorfun,colorminv=colorminv,colormaxv=colormaxv)

"Plot `Mean` and `SEM` of responses for each condition"
function plotcondresponse(rs,ctc;factors=propertynames(ctc),u=0,color=:auto,style=:path,title="",projection=:none,grid=false,linewidth=:auto,legend=:best,response=[])
    plotcondresponse(Dict(u=>rs),ctc,factors;color,style,title,projection,linewidth,legend,response,grid)
end
function plotcondresponse(urs::Dict,ctc::DataFrame,factors;color=huecolors(length(urs)),style=:path,grid=false,projection=:none,title="",linewidth=:auto,legend=:best,response=[])
    df = condresponse(urs,ctc,factors)
    plotcondresponse(df;color,style,title,projection,linewidth,legend,response,grid)
end
function plotcondresponse(urs::Dict,cond::DataFrame;color=huecolors(length(urs)),style=:path,projection=:none,title="",grid=false,linewidth=:auto,legend=:best,response=[])
    df = condresponse(urs,cond)
    plotcondresponse(df;color,style,title,projection,linewidth,legend,response,grid)
end
function plotcondresponse(mseuugc::DataFrame;group=:ug,color=:auto,style=:path,projection=:none,title="",grid=false,
                            linestyle=:solid,linewidth=:auto,legend=:best,response=[],responsetype=:Response)
    # sort to the group order plot uses
    gs = unique(mseuugc[!,[:u,:ug]])
    if size(gs,1)>1
        gi = sortperm(gs);gs=gs[gi,:]
        color isa Vector && (color=permutedims(color[gi]))
        linewidth isa Vector && (linewidth=permutedims(linewidth[gi]))
        linestyle isa Vector && (linestyle=permutedims(linestyle[gi]))
    end
    factors=setdiff(propertynames(mseuugc),[:m,:se,:u,:ug])
    nfactor=length(factors)
    if nfactor==1
        factor=factors[1]
        # if typeof(mseuc[!,factor][1]) <: Array
        #     map!(string,mseuc[factor],mseuc[factor])
        #     style=:bar
        # end
    elseif nfactor==2
        fm,fse,fa = factorresponse(mseuugc)
        maxlim = skipmissing(fm) |> extrema .|> abs |> maximum
        minlim = skipmissing(fm) |> minimum
        yfactor,xfactor = collect(keys(fa))
        y,x = collect(values(fa))
        clims = minlim < 0 ? (-maxlim,maxlim) : (0,maxlim)
    else
        mseuugc[:Condition]=condstring(mseuugc[!,factors])
        factor=:Condition
        style=:bar
    end
    if nfactor==2
        heatmap(x,y,fm;color,title,legend,xaxis=(factorunit(xfactor)),yaxis=(factorunit(yfactor)),colorbar_title=factorunit(responsetype),clims)
    else
        if projection==:polar # close curve
            c0 = mseuugc[mseuugc[!,factor].==0,:]
            c0[!,factor].=360
            mseuugc = [mseuugc;c0]
            mseuugc[!,factor]=deg2rad.(mseuugc[!,factor])
            be = backend_name()
            if be == :gr
                legend=(0.86,0.96)
            elseif be == :pyplot
                legend=(0.91,0.81)
            end
        end
        sort!(mseuugc,factor)
        if projection==:polar
            if be == :gr
                p = @df mseuugc plot(cols(factor),:m;yerror=:se,group=cols(group),line=style,markerstrokecolor=color,color,
                label=permutedims(["$(k.ug)$(k.u)" for k in eachrow(gs)]),grid,projection,legend,
                xaxis=(factorunit(factor)),yaxis=(factorunit(responsetype)),title,linestyle,linewidth)
            elseif be == :pyplot
                p = @df mseuugc plot(cols(factor),:m;ribbon=:se,group=cols(group),line=style,markerstrokecolor=color,color,
                label=permutedims(["$(k.ug)$(k.u)" for k in eachrow(gs)]),grid,projection,legend,
                xticks=0:0.25π:1.75π,xformatter=x->round(Int,rad2deg(x)),title,linestyle,linewidth)
            end
        else
            p = @df mseuugc plot(cols(factor),:m;ribbon=:se,group=cols(group),line=style,markerstrokecolor=color,color,
            label=permutedims(["$(k.ug)$(k.u)" for k in eachrow(gs)]),grid,projection,legend,
            xaxis=(factorunit(factor)),yaxis=(factorunit(responsetype)),title,linestyle,linewidth)
        end
        if !isempty(response)
            for i in response
                hline!(p,[i[1]],ribbon=[i[2]],color=colors,legend=false)
            end
        end
        p
    end
end

plotpsth(rvv::RVVector,binedges::RealVector;timeline=[0],colors=[:auto],title="")=plotpsth(rvv,binedges,DataFrame(Factor="Value",i=[1:length(rvv)]),timeline=timeline,colors=colors,title=title)
function plotpsth(rvv::RVVector,binedges::RealVector,ctc::DataFrame,factor;timeline=[0],colors=nothing,title="")
    msexc = psth(rvv,binedges,ctc,factor)
    plotpsth(msexc,timeline=timeline,colors=colors==nothing ? huecolors(length(levels(msexc[:c]))) : colors,title=title)
end
function plotpsth(rvv::RVVector,binedges::RealVector,cond::DataFrame;timeline=[0],colors=huecolors(nrow(cond)),title="")
    msexc = psth(rvv,binedges,cond)
    plotpsth(msexc,timeline=timeline,colors=colors,title=title)
end
function plotpsth(msexc::DataFrame;timeline=[0],colors=[:auto],title="")
    @df msexc Plots.plot(:x,:m,ribbon=:se,group=:c,fillalpha=0.2,color=reshape(colors,1,:))
    vline!(timeline,line=(:grey),label="TimeLine",grid=false,xaxis=(factorunit(:Time)),yaxis=(factorunit(responsetype)),title=(title))
end

function plotsta(ps;sizepx=size(ps),sizedeg=nothing,ppd=45,index=nothing,filter=Kernel.gaussian(1),title="",color="redblue",r=[extrema(ps)...],bg="white")
    nd = ndims(ps)
    if nd==1
        if isnothing(index)
            sta = reshape(ps,sizepx)
        else
            sta = fill(mean(ps),sizepx)
            sta[index] = ps
        end
    elseif nd==2
        sta = ps
    end
    if !isnothing(filter)
        sta = imfilter(sta,filter)
    end
    if !isnothing(sizedeg)
        ppd = first(sizepx./sizedeg)
    end

    x = vec([(j-1)/ppd for i in 1:sizepx[1], j in 1:sizepx[2]])
    y = vec([(i-1)/ppd for i in 1:sizepx[1], j in 1:sizepx[2]])
    z = vec(sta)
    if !isnothing(index)
        x=x[index];y=y[index];z=z[index]
    end
    xlim=[extrema(x)...];ylim=[extrema(y)...];psize=2;width=psize*length(unique(x));height=psize*length(unique(y))

    DataFrame(z=z,x=x,y=y) |> @vlplot(mark={:rect,size=psize,strokeWidth=0},width=width,height=height,background=bg,
    x={"x:o",title="X (deg)",axis={values=xlim,format=".1",labelAngle=0}},
    y={"y:o",title="Y (deg)",axis={values=ylim,format=".1"}},
    color={"z:q",title="",scale={domain=r}},title={text=title},
    config={range={heatmap={scheme=color,extent=[1,0]}},
    view={strokeWidth=0,fill=bg}})

    # x=[(i-1)/ppd for i in 1:imagesize[1]]
    #
    # color=:coolwarm
    # if true
    #     cg = cgrad(color)
    #     minv,maxv = extrema(sta)
    #     vr=maxv-minv
    #     sta = map(i->coloralpha(RGB(cg[(i-minv)/vr]),0),sta)
    #     sta[index] = coloralpha.(sta[index],1)
    # end
    # Plots.heatmap(x,x,sta,color=:coolwarm,ratio=:equal,yflip=true,leg=true,framestyle=:grid,title=title,
    # xlabel="Position_X (Deg)",ylabel="Position_Y (Deg)",xtick=xlim,ytick=[])
end

"Plot Analog Signals"
function plotanalog(data;x=nothing,y=nothing,fs=0,xext=0,timeline=[],xlabel="Time",ylabel="Depth",clims=nothing,
    xunit=:ms,yunit=:μm,cunit=:v,plottype=:heatmap,hy=0,color=:coolwarm,layer=nothing,n=nothing,aspectratio=:auto)
    nd=ndims(data);ds=1
    if nd==1
        if isnothing(x)
            x=1:length(data)
            if fs>0
                x = x./fs.-xext
                xunit==:ms && (x*=1000)
            end
        end
        if isnothing(clims)
            if cunit==:v
                lim = maximum(abs.(data))
                clims = (-lim,lim)
            elseif cunit==:uv
                ds = 1e6
                lim = maximum(abs.(data))*ds
                clims = (-lim,lim)
            elseif cunit == :fr
                lim = maximum(data)
                clims = (0,lim)
            else
                clims = :auto
            end
        end
        xlabel = factorunitstring(xlabel,xunit)
        ylabel = factorunitstring(ylabel,yunit)
        p=plot(x,data*ds;legend=false,color_palette=color,grid=false,ylims=clims,xlabel,ylabel,aspectratio)
    elseif nd==2
        if isnothing(x)
            x=1:size(data,2)
            if fs>0
                x = x./fs.-xext
                xunit==:ms && (x*=1000)
            end
        end
        if isnothing(clims)
            if cunit==:v
                lim = maximum(abs.(data))
                clims = (-lim,lim)
            elseif cunit==:uv
                ds = 1e6
                lim = maximum(abs.(data))*ds
                clims = (-lim,lim)
            elseif cunit == :fr
                lim = maximum(data)
                clims = (0,lim)
            else
                clims = :auto
            end
        end
        xlabel = factorunitstring(xlabel,xunit)
        ylabel = factorunitstring(ylabel,yunit)
        if plottype==:heatmap
            if isnothing(y)
                y = (1:size(data,1))
                hy>0 && (y*=hy)
            end
            p=heatmap(x,y,data*ds;color=color,clims=clims,xlabel,ylabel,aspectratio,xlim=extrema(x),ylim=extrema(y))
            if !isnothing(n)
                pn = clampscale(n) .* (maximum(x)-minimum(x)) .* 0.2 .+ minimum(x)
                plot!(p,pn,y,color=:seagreen,leg=false,xlim=extrema(x),ylim=extrema(y))
            end
        else
            p=plot(x,data'*ds;legend=false,color_palette=color,grid=false,ylims=clims,xlabel,ylabel,aspectratio)
        end
    end
    isempty(timeline) || vline!(p,timeline,line=(:grey),label="TimeLine",leg=false)
    if !isnothing(layer)
        anno = [(minimum(x)+3,mean(layer[k]),text(k,6,:gray10,:center,:left)) for k in keys(layer)]
        hline!(p,[l[2] for l in values(layer)],linecolor=:gray25,legend=false,lw=0.5,annotations=anno)
    end
    return p
end

plotunitposition(spike::Dict;layer=nothing,color=nothing,alpha=0.4,title="",markersize=5,unitidsize=3,size=(600,450)) = plotunitposition(spike["unitposition"];unitgood=spike["unitgood"],chposition=spike["chposition"],unitid=spike["unitid"],layer,color,alpha,title,markersize,unitidsize,size)
function plotunitposition(unitposition;unitgood=[],chposition=[],unitid=[],layer=nothing,color=nothing,alpha=0.4,title="",markersize=5,unitidsize=3,size=(600,450))
    nunit = Base.size(unitposition,1);ngoodunit = isempty(unitgood) ? nunit : count(unitgood);ustr = "$ngoodunit / $nunit"
    xlims = isempty(chposition) ? (minimum(unitposition[:,1])-4,maximum(unitposition[:,1])+2) : (minimum(chposition[:,1])-5,maximum(chposition[:,1])+5)
    p = plot(;legend=:topright,xlabel="X (μm)",ylabel="Y (μm)",xlims,size)
    if !isempty(chposition)
        scatter!(p,chposition[:,1],chposition[:,2],markershape=:rect,markerstrokewidth=0,markersize=1.5,color=:grey70,label="Electrode")
    end
    if isnothing(color)
        if isempty(unitgood)
            color = :gray30
        else
            color = map(i->i ? :darkgreen : :gray30,unitgood)
        end
    end
    if !isnothing(layer)
        lx = xlim[1]+2
        hline!(p,[layer[k][1] for k in keys(layer)],linestyle=:dash,annotations=[(lx,layer[k][1],text(k,5,:gray20,:bottom)) for k in keys(layer)],linecolor=:gray30,legend=false)
    end
    if isempty(unitid)
        scatter!(p,unitposition[:,1],unitposition[:,2];label=ustr,color,alpha,markerstrokewidth=0,markersize,title)
    else
        scatter!(p,unitposition[:,1],unitposition[:,2];label=ustr,color,alpha,markerstrokewidth=0,markersize,series_annotations=text.(unitid,unitidsize,:gray10,:center),title)
    end
    return p
end

function plotunitpositionproperty(unitposition;ori=nothing,os=nothing,dir=nothing,ds=nothing,sf=nothing,width=500,height=400,title="",layer=nothing)
    df = DataFrame(x=unitposition[:,1],y=unitposition[:,2],m=Any[:circle for _ in 1:size(unitposition,1)],a=0.0,sw=0.0,c=0.0,s=10.0)
    l = DataFrame()
    xlim = [minimum(df[:,:x])-4,maximum(df[:,:x])+2]
    ylim = [minimum(df[:,:y])-100,maximum(df[:,:y])+100]
    if !isnothing(sf)
        df[!,:c]=sf
    end
    if !isnothing(ori)
        df[!,:m] .= :stroke
        df[!,:a] = -ori
        df[!,:sw] .= 1
        if !isnothing(os)
            df[!,:s]=os
        else
            df[!,:s].=160
        end
    end
    if !isnothing(dir)
        t = copy(df)
        arrowpath = "M 0 -0.1 H 1 V -0.3 L 1.6 0 L 1 0.3 V 0.1 H 0 Z"
        t[!,:m] .= arrowpath
        t[!,:a] = -dir
        t[!,:sw] .= 0
        if !isnothing(ds)
            t[!,:s]=ds
        else
            t[!,:s].=200
        end
        df = [df;t]
    end
    if !isnothing(layer)
        l[!,:x] = fill(xlim[1],length(layer))
        l[!,:x2] .= xlim[2]
        l[!,:y] = [v[1] for v in values(layer)]
        l[!,:y2] = l[!,:y]
        l[!,:l] = collect(keys(layer))
        ylim = [extrema([ylim;l[:,:y]])...].+[-50,100]
    end
    @vgplot(height=height,width=width,padding=5,data=[:df=>df,:l=>l],
    marks=[
    {
    type="rule",
    from={data="l"},
    encode={
        update={
        x={field="x",scale="x"},
        y={field="y",scale="y"},
        x2={field="x2",scale="x"},
        y2={field="y2",scale="y"},
        strokeDash={value=[4,2]},
        strokeWidth={value=0.5},
        stroke={value="dimgray"}
        }}
    },
    {
    type="text",
    from={data="l"},
    encode={
        update={
        x={field="x",scale="x"},
        y={field="y",scale="y"},
        text={field="l"},
        align={value="left"},
        baseline={value="bottom"},
        fontSize={value=7},
        dx={value=15}
        }}
    },
    {
    type="symbol",
    from={data="df"},
    encode={
        update={
        x={field="x",scale="x"},
        y={field="y",scale="y"},
        shape={field="m"},
        angle={field="a"},
        size={field="s",scale="s"},
        strokeWidth={field="sw"},
        stroke={field="c",scale="c"},
        fill={field="c",scale="c"}
        }}
    }
    ],
    scales=[
    {
        name="x",
        nice=false,
        zero=false,
        range="width",
        domain=xlim,
        type="linear",
        round=true
    },
    {
        name="y",
        nice=false,
        zero=false,
        range="height",
        domain=ylim,
        type="linear",
        round=true
    },
    {
        name="c",
        nice=false,
        zero=false,
        round=false,
        type="linear",
        range={scheme = "plasma",extent=[0.8,0.2]},
        domain={data="df",field="c"}
    },
    {
        name="s",
        nice=false,
        zero=false,
        domain={data="df",field="s"},
        type="linear",
        round=true,
        range=[10, 200]
    }
    ],
    axes=[
    {
        domain=true,
        tickCount=5,
        grid=false,
        title="Position_X (μm)",
        scale="x",
        orient="bottom"
    },
    {
        domain=true,
        tickCount=5,
        grid=false,
        title="Position_Y (μm)",
        scale="y",
        orient="left"
    }
    ],
    title={
    text = title
    },
    legends=[
    {
    type="gradient",
    fill="c",
    title="SF"
    },
    {
    type="symbol",
    symbolType="stroke",
    size="s",
    title="1-CV"
    }
    ])
end

function plotunitpositionimage(unitposition,image;width=800,height=600,markersize=20,title="",layer=nothing)
    tempimagedir = "tempimage";nu=size(unitposition,1)
    isdir(tempimagedir) || mkpath(tempimagedir)
    urls = map(i->"tempimage/$(i).png",1:nu)
    foreach(i->save(urls[i],image[i]),1:nu)

    df = DataFrame(x=unitposition[:,1],y=unitposition[:,2],m=urls,s=markersize)
    l = DataFrame()
    # xlim = [minimum(df[:,:x])-4,maximum(df[:,:x])+2]
    # ylim = [minimum(df[:,:y])-100,maximum(df[:,:y])+100]
    xlim = [minimum(df[:,:x])-0.2,maximum(df[:,:x])+0.2]
    ylim = [minimum(df[:,:y])-0.2,maximum(df[:,:y])+0.2]
    # xlim = [minimum(df[:,:x])-2,maximum(df[:,:x])+2]
    # ylim = [minimum(df[:,:y])-2,maximum(df[:,:y])+2]
    if !isnothing(layer)
        l[!,:x] = fill(xlim[1],length(layer))
        l[!,:x2] .= xlim[2]
        l[!,:y] = [v[1] for v in values(layer)]
        l[!,:y2] = l[!,:y]
        l[!,:l] = collect(keys(layer))
        ylim = [extrema([ylim;l[:,:y]])...].+[-50,100]
    end
    @vgplot(height=height,width=width,padding=5,data=[:df=>df,:l=>l],
    marks=[
    {
    type="rule",
    from={data="l"},
    encode={
        update={
        x={field="x",scale="x"},
        y={field="y",scale="y"},
        x2={field="x2",scale="x"},
        y2={field="y2",scale="y"},
        strokeDash={value=[4,2]},
        strokeWidth={value=0.5},
        stroke={value="dimgray"}
        }}
    },
    {
    type="text",
    from={data="l"},
    encode={
        update={
        x={field="x",scale="x"},
        y={field="y",scale="y"},
        text={field="l"},
        align={value="left"},
        baseline={value="bottom"},
        fontSize={value=7},
        dx={value=15}
        }}
    },
    {
    type="image",
    from={data="df"},
    encode={
        update={
        x={field="x",scale="x"},
        y={field="y",scale="y"},
        url={field="m"},
        width={field="s"},
        height={field="s"},
        align={value="center"},
        baseline={value="middle"}
        }}
    }
    ],
    scales=[
    {
        name="x",
        nice=false,
        zero=false,
        range="width",
        domain=xlim,
        type="linear",
        round=true
    },
    {
        name="y",
        nice=false,
        zero=false,
        range="height",
        domain=ylim,
        type="linear",
        round=true
    }
    ],
    axes=[
    {
        domain=false,
        tickCount=0,
        grid=false,
        # title="Position_X (μm)",
        title="",
        scale="x",
        orient="bottom"
    },
    {
        domain=false,
        tickCount=0,
        grid=false,
        # title="Position_Y (μm)",
        title="",
        scale="y",
        orient="left"
    }
    ],
    title={
    text = title
    }
    )
end

function plotunitlayerimage(unitlayer,image;width=800,height=600,markersize=40,title="",unitid=nothing)
    tempimagedir = joinpath(pwd(),"tempimage");nu=length(unitlayer)
    isdir(tempimagedir) || mkpath(tempimagedir)
    urls = map(i->"tempimage/$(i).png",1:nu)
    foreach(i->save(urls[i],image[i]),1:nu)

    x = zeros(nu)
    for r in eachrow(condin(DataFrame(l=unitlayer)))
        x[r.i] = 1:r.n
    end
    df = DataFrame(x=x,y=unitlayer,m=urls)
    if isnothing(unitid)
        df[:u] = ""
    else
        df[:u] = unitid
    end

    @vlplot(height=height,width=width,padding=5,title=title,data=df,
    y={"y:o",axis={title="Layer"}},
    x={"x",axis={grid=false,title=nothing,ticks=false,domain=false,labels=false}}) +
    @vlplot(mark={:image,width=markersize,height=markersize,align=:center,baseline=:middle}, url=:m) +
    @vlplot(mark={:text,align=:center,baseline=:bottom,dy=30,fontSize=7},text=:u)
end

function plotcircuit(unitposition,unitid,projs;unitgood=[],eunits=[],iunits=[],projweights=[],layer=nothing,showuid=true,showmode=:none)
    vpi = indexin(unique(projs),projs)
    projs=projs[vpi];np=length(projs)
    if !isempty(projweights)
        projweights=projweights[vpi]
        t = abs.(projweights)
        pcolor = coloralpha.(RGB(0.5,0.5,0.5),t./maximum(t))
    else
        pcolor=:gray50
    end
    nu = length(unitid)
    nsu = isempty(unitgood) ? nu : count(unitgood)
    uidi = Dict(i=>findfirst(i.==unitid) for i in unitid)
    nsupair=binomial(nsu,2);ss = "$nsu: $np/$nsupair($(round(np/nsupair*100,digits=3))%)"
    xlim = (minimum(unitposition[:,1])-5,maximum(unitposition[:,1])+5)
    ylim = (minimum(unitposition[:,2]),maximum(unitposition[:,2]))
    p = plot(legend=:topright,xlabel="Position_X (um)",ylabel="Position_Y (um)",xlims=xlim,grid=false)

    for i in projs
        t=hcat(unitposition[uidi[i[1]],:],unitposition[uidi[i[2]],:])
        plot!(p,t[1,:],t[2,:],linewidth=0.3,color=pcolor,arrow=arrow(:closed,:head,0.45,0.12))
    end

    color = fill(:gray30,nu)
    color[unitgood] .= :darkgreen
    if !isempty(eunits)
        color[map(i->uidi[i],eunits)] .= :darkred
    end
    if !isempty(iunits)
        color[map(i->uidi[i],iunits)] .= :darkblue
    end
    if showmode == :su
        showuidi = unitgood
    elseif showmode == :circuit
        cuid=[]
        foreach(p->append!(cuid,[p...]),projs)
        showuidi=indexin(unique(cuid),unitid)
    else
        showuidi = 1:nu
    end
    if showuid
        scatter!(p,unitposition[showuidi,1],unitposition[showuidi,2],label=ss,color=color[showuidi],alpha=0.4,markerstrokewidth=0,markersize=6,series_annotations=text.(unitid[showuidi],3,:gray10,:center),legend=false)
    else
        scatter!(p,unitposition[showuidi,1],unitposition[showuidi,2],label=ss,color=color[showuidi],alpha=0.4,markerstrokewidth=0,markersize=5,legend=false)
    end
    if !isnothing(layer)
        hline!(p,[layer[k][1] for k in keys(layer)],linestyle=:dash,annotations=[(xlim[1]+2,layer[k][1],text(k,5,:gray20,:bottom)) for k in keys(layer)],linecolor=:gray30,legend=false)
    end
    annotate!(p,[(xlim[2]-6,ylim[2]-100,text(ss,6,:gray20,:bottom))])
    return p
end

"plot hartley subspace"
function plothartleysubspace(ps,nk,dk;color=:grays)
    n = 2nk+1
    p=plot(layout=(n,n),size=(120n,120n),leg=false,clims=(-1,1),frame=:none,aspect_ratio=:equal)
    xy=0:0.01:1
    for (kx,ky,phase) in ps
        r = -Int(ky/dk)+nk+1; c = Int(kx/dk)+nk+1
        cg = [cas(i,j,kx=kx,ky=ky,phase=phase) for j in xy,i in xy]
        heatmap!(p,subplot=c+n*(r-1),cg,color=color)
    end
    p
end

# function savefig(fig,filename::AbstractString;path::AbstractString="",format::AbstractString="svg")
#     f = joinpath(path,"$filename.$format")
#     if !ispath(path)
#         mkpath(path)
#     end
#     if format=="svg"
#         format = "$format+xml"
#     end
#     open(f, "w") do io
#         writemime(io,"image/$format",fig)
#     end
# end
#
# function savefig(fig::Gadfly.Plot,filename::AbstractString;path::AbstractString="",format::AbstractString="svg",width=22cm,height=13cm,dpi=300)
#     f = joinpath(path,"$filename.$format")
#     if !ispath(path)
#         mkpath(path)
#     end
#     if format=="svg"
#         format = SVG(f,width,height)
#     elseif format=="png"
#         format = PNG(f,width,height,dpi=dpi)
#     elseif format=="pdf"
#         format = PDF(f,width,height,dpi=dpi)
#     end
#     draw(format,fig)
# end

ext: jl


text:
using NLsolve, StaticArrays

function f!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function j!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

nlsolve(f!, j!, @MArray [ 0.1; 1.2], )


text:
"""
    FillMap(λ, (m, n))::FillMap
    FillMap(λ, m, n)::FillMap

Construct a (lazy) representation of an operator whose matrix representation
would be an m×n-matrix filled constantly with the value `λ`.
"""
struct FillMap{T} <: LinearMap{T}
    λ::T
    size::Dims{2}
    function FillMap(λ::T, dims::Dims{2}) where {T}
        (dims[1]>=0 && dims[2]>=0) || throw(ArgumentError("dims of FillMap must be non-negative"))
        return new{T}(λ, dims)
    end
end

FillMap(λ, m::Int, n::Int) = FillMap(λ, (m, n))

# properties
Base.size(A::FillMap) = A.size
MulStyle(A::FillMap) = FiveArg()
LinearAlgebra.issymmetric(A::FillMap) = A.size[1] == A.size[2]
LinearAlgebra.ishermitian(A::FillMap) = isreal(A.λ) && A.size[1] == A.size[2]
LinearAlgebra.isposdef(A::FillMap) = (size(A, 1) == size(A, 2) == 1 && isposdef(A.λ))
Base.:(==)(A::FillMap, B::FillMap) = A.λ == B.λ && A.size == B.size

LinearAlgebra.adjoint(A::FillMap) = FillMap(adjoint(A.λ), reverse(A.size))
LinearAlgebra.transpose(A::FillMap) = FillMap(transpose(A.λ), reverse(A.size))

function Base.:(*)(A::FillMap, x::AbstractVector)
    T = typeof(oneunit(eltype(A)) * (zero(eltype(x)) + zero(eltype(x))))
    return fill(iszero(A.λ) ? zero(T) : A.λ*sum(x), A.size[1])
end

function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector)
    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
end

_unsafe_mul!(Y::AbstractMatrix, A::FillMap, x::Number) = fill!(Y, A.λ*x)
function _unsafe_mul!(Y::AbstractMatrix, A::FillMap, x::Number, α::Number, β::Number)
    LinearAlgebra._rmul_or_fill!(Y, β)
    Y .+= A.λ*x*α
    return Y
end

function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector, α::Number, β::Number)
    if iszero(α)
        !isone(β) && rmul!(y, β)
    else
        temp = A.λ * sum(x) * α
        if iszero(β)
            y .= temp
        elseif isone(β)
            y .+= temp
        else
            y .= y .* β .+ temp
        end
    end
    return y
end

Base.:(+)(A::FillMap, B::FillMap) = A.size == B.size ? FillMap(A.λ + B.λ, A.size) : throw(DimensionMismatch())
Base.:(-)(A::FillMap) = FillMap(-A.λ, A.size)
Base.:(*)(λ::Number, A::FillMap) = FillMap(λ * A.λ, size(A))
Base.:(*)(A::FillMap, λ::Number) = FillMap(A.λ * λ, size(A))
Base.:(*)(λ::RealOrComplex, A::FillMap) = FillMap(λ * A.λ, size(A))
Base.:(*)(A::FillMap, λ::RealOrComplex) = FillMap(A.λ * λ, size(A))

function Base.:(*)(A::FillMap, B::FillMap)
    check_dim_mul(A, B)
    return FillMap(A.λ*B.λ*size(A, 2), (size(A, 1), size(B, 2)))
end


text:
using DashBootstrapComponents

row = html_div([
    dbc_row(dbc_col(html_div("A single column"))),
    dbc_row([
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns")),
    ]),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: Hairpin.mso  Z bicca $
#*------------------------------------------------------------------
type Hairpin_LMTD
	Hairpin_LMTD()=begin
		new(
			Hairpin_Basic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flow Direction",
				:Valid=>["counter","cocurrent"],
				:Default=>"cocurrent"
			)),
			LMTD_Basic(Dict{Symbol,Any}(
				:Brief=>"LMTD Method of Calculation",
				:Symbol=>" "
			)),
			[
				:(_base_1.Details.Q = _base_1.Details.Ud*_base_1.Pi*_base_1.DoInner*(2*_base_1.Lpipe)*Method.LMTD),
				:(_base_1.Details.Q = _base_1.Qestimated),
				:(Method.Fc = 1),
				:(Method.DT0 = _base_1.InletOuter.T - _base_1.InletInner.T),
				:(Method.DTL = _base_1.OutletOuter.T - _base_1.OutletInner.T),
				:(Method.DT0 = _base_1.InletOuter.T - _base_1.OutletInner.T),
				:(Method.DTL = _base_1.OutletOuter.T - _base_1.InletInner.T),
				:(Method.DT0 = _base_1.InletInner.T - _base_1.InletOuter.T),
				:(Method.DTL = _base_1.OutletInner.T - _base_1.OutletOuter.T),
				:(Method.DT0 = _base_1.InletInner.T - _base_1.OutletOuter.T),
				:(Method.DTL = _base_1.OutletInner.T - _base_1.InletOuter.T),
			],
			[
				"Duty","Duty  Estimated","LMTD Correction Factor - True counter ou cocurrent flow","Temperature Difference at Inlet - Cocurrent Flow","Temperature Difference at Outlet - Cocurrent Flow","Temperature Difference at Inlet - Counter Flow","Temperature Difference at Outlet - Counter Flow","Temperature Difference at Inlet - Cocurrent Flow","Temperature Difference at Outlet - Cocurrent Flow","Temperature Difference at Inlet - Counter Flow","Temperature Difference at Outlet - Counter Flow",
			],
			[:FlowDirection,],
			[:Method,]
		)
	end
	_base_1::Hairpin_Basic
	FlowDirection::DanaSwitcher
	Method::LMTD_Basic
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Hairpin_LMTD
function setEquationFlow(in::Hairpin_LMTD)
	let switch=CalculationApproach
		if switch=="Full"
			addEquation(1)
		elseif switch=="Simplified"
			addEquation(2)
		end
	end
	addEquation(3)
	let switch=HotSide
		if HotSide==InletInner.T > InletOuter.T
			set(switch,"inner")
		end
		if HotSide==InletInner.T < InletOuter.T
			set(switch,"outer")
		end
		if switch=="outer"
			let switch=FlowDirection
				if switch=="cocurrent"
					addEquation(4)
					addEquation(5)
				elseif switch=="counter"
					addEquation(6)
					addEquation(7)
				end
			end
		elseif switch=="inner"
			let switch=FlowDirection
				if switch=="cocurrent"
					addEquation(8)
					addEquation(9)
				elseif switch=="counter"
					addEquation(10)
					addEquation(11)
				end
			end
		end
	end
end
function atributes(in::Hairpin_LMTD,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Icon]="icon/hairpin"
	fields[:Pallete]=true
	fields[:Brief]="Hairpin Heat Exchanger - LMTD Method"
	fields[:Info]="to be documented."
	drive!(fields,_)
	return fields
end
Hairpin_LMTD(_::Dict{Symbol,Any})=begin
	newModel=Hairpin_LMTD()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Example of composite type
mutable struct Person
    first_name
    last_name
    age::Integer
 end

 # We create an instance of Person
 person = Person("John","Doe",25)
 print(person)

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
# This file uses MathML.jl and libsbml.jl to convert SBML models to ReactionSystems

struct SbmlModel
    # Basic Functionality
    parameters#::Array{Pair{Variable,Float64}}  # [parameter=>value, ...]
    compartments#::Array{Pair{Variable,Float64}}  # [nucleus=>1.0,cytoplasm=>2.0,...]
    species#::Array{Pair{Variable,Tuple{Float64,Variable}}}  # [specie=>(u0,compartment), ...]
    reactions#::Array{tuple{SymbolicUtils.Add,Array{Variable},Array{Variable},Array{Int},Array{Int}}}
         # [(fullkineticlaw,[subtrate,...],[product,...],[substoich,...],[prodstoich,...]),...]
    ## Future components
    # maybe add fields to store info about events, piecewise simulation, conversionFactors etc.
end
SbmlModel(file::String) = SbmlModel(readxml(file))
SbmlModel(doc::EzXML.Document) = _process_doc(doc)

create_var(x) = Num(Variable(Symbol(x)))
# create_var(x, iv) = Num(Sym{FnType{Tuple{Real}}}(Symbol(x))(Variable(Symbol(iv)))).val
# create_var(x, iv) = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(x)))(Variable(Symbol(iv)))
create_param(x) = Num(Sym{ModelingToolkit.Parameter{Real}}(Symbol(x)))

""" Convert SBML document into SbmlModel """
function _process_doc(doc)
    # doc = make_extensive(doc)
    doc = promotelocalparameters(doc)
    sysinfo = sbml_to_sysinfo(doc)
    pars = vcat(sysinfo["listOfParameters"],
                sysinfo["listOfReactions/x:reaction/x:kineticLaw/x:listOfParameters"],
                sysinfo["listOfReactions/x:reaction/x:kineticLaw/x:listOfLocalParameters"])
    pars = build_par_map(pars)
    comps = build_comp_map(sysinfo["listOfCompartments"])
    spec = build_spec_map(sysinfo["listOfSpecies"])
    substitutions = _mathml_substitutions(sysinfo)
    reactions = build_reactions(sysinfo["listOfReactions"], substitutions)
    SbmlModel(pars,comps,spec,reactions)
end

function _mathml_substitutions(sysinfo)
    pars = vcat(sysinfo["listOfCompartments"],
                sysinfo["listOfParameters"],
                sysinfo["listOfReactions/x:reaction/x:kineticLaw/x:listOfParameters"],
                sysinfo["listOfReactions/x:reaction/x:kineticLaw/x:listOfLocalParameters"])
    vars = sysinfo["listOfSpecies"]
    par_subs = _par_subs(pars)
    # var_subs = _var_subs(vars)
    Dict(par_subs)
end

function _par_subs(nodes)
    mathml = @. create_var(getindex(nodes, "id"))
    sbmlmodel = @. create_param(getindex(nodes, "id"))
    mathml .=> sbmlmodel
end

#=function _var_subs(nodes)
    mathml = @. Num(Variable(Symbol(getindex(nodes, "id"))))
    sbmlmodel = @. Num(Sym{ModelingToolkit.Parameter{Real}}(Symbol(getindex(nodes, "id"))))
    mathml .=> sbmlmodel
end=#

""" Extract kineticLaws, Reactants, Products and their stoichiometry from SBML """
function build_reactions(listofreactions::Vector{EzXML.Node}, substitutions)
    reactions = Tuple{Num,Array{Num,1},Array{Num,1},Array{Int64,1},Array{Int64,1}}[]
    for reaction in listofreactions
        reactants, r_stoich = _getlistofspeciesreference(reaction,"Reactant")
        products, p_stoich = _getlistofspeciesreference(reaction,"Product")
        kineticlaw = getkineticlaw(reaction)
        kineticlaw = parse_node(getmath(kineticlaw))[1]
        kineticlaw = substitute(kineticlaw, substitutions)
        thisreaction = [Tuple{Num,Array{Num,1},Array{Num,1},Array{Int64,1},Array{Int64,1}}((kineticlaw,reactants,products,r_stoich,p_stoich))]
        append!(reactions,thisreaction)
    end
    reactions
end

""" Extract species and their stroichiometry from an SBML reaction """
function _getlistofspeciesreference(reaction::EzXML.Node,type)
    if nodename(reaction) != "reaction"
        @error("Input node must be a reaction but is $(nodename(reaction)).")
    end
    listnodes = [node for node in eachelement(reaction) if nodename(node) == "listOf$(type)s"]
    if length(listnodes) > 1
        @error("SBML files with reactions with more than one listOf$(type)s are not supported.")
    end
    if length(listnodes) == 0
        return (nothing, nothing)
    end
    listnode = listnodes[1]
    spec = [create_var(getindex(node, "species")) for node in eachelement(listnode) if nodename(node) == "speciesReference"]
    stoich = [Int(Meta.parse(getindex(node, "stoichiometry"))) for node in eachelement(listnode) if nodename(node) == "speciesReference"]
    (spec, stoich)
end

""" Extract kineticLaw node from an SBML reaction """
function getkineticlaw(reaction::EzXML.Node)
    kineticlaws = [node for node in eachelement(reaction) if nodename(node) == "kineticLaw"]
    if length(kineticlaws) > 1
        @error("Reactions with more than one kinticLaw are not supported.")
    end
    kineticlaw = kineticlaws[1]
end

""" Extract math node from an SBML kineticLaw """
function getmath(kineticlaw::EzXML.Node)
    maths = [node for node in eachelement(kineticlaw) if nodename(node) == "math"]
    if length(maths) > 1
        @error("kineticLaws with more than one math node are not supported.")
    end
    math = maths[1]
end

""" Extract parameters from SBML """
function build_par_map(parnodes::Vector{EzXML.Node})
    ids = @. create_param(getindex(parnodes, "id"))
    vals = @. Float64(Meta.parse(getindex(parnodes, "value")))
    ids .=> vals
end

""" Extract compartments from SBML """
function build_comp_map(compnodes::Vector{EzXML.Node})
    ids = @. create_param(getindex(compnodes, "id"))
    sizes = @. Float64(Meta.parse(getindex(compnodes, "size")))
    ids .=> sizes
end

""" Extract species from SBML """
function build_spec_map(specnodes::Vector{EzXML.Node})
    # ids = @. Num(Variable{Float64}(Symbol(getindex(specnodes, "id"))))
    ids = @. create_var(getindex(specnodes, "id"))
    inits = @. Float64(Meta.parse(getindex(specnodes, "initialAmount")))
    comps = @. create_param(getindex(specnodes, "compartment"))
    ids .=> tuple.(inits,comps)
end

#=function expand_reversible(doc::EzXML.Node)
    doc = deepcopy(doc)
    root = doc.root
    ns = namespace(root)
    reactions = findall("//x:reaction", root, ["x"=>ns])
    for reaction in reactions
        kineticlaws = [node for node in eachelement(reaction) if nodename(node) == "kineticLaw"]
        if length(kineticlaws) > 1
            @error("SBML reactions with more than one kinetic law are not supported.")
        end
        kineticlaw = kineticlaws[1]
        locallistsofparameters = [node for node in eachelement(kineticlaw) if occursin("Parameters",nodename(node))]
        locallistofparameters = [node for list in locallistsofparameters for node in eachelement(list)]
        for par in locallistofparameters
            par["id"] = reaction["id"]*"_"*par["id"]
            par["name"] = reaction["name"]*"_"*par["name"]
        end
    end
    doc
end=#

""" Assign gloablly unique ids to local parameters """
function promotelocalparameters(doc::EzXML.Document)
    doc = deepcopy(doc)
    root = doc.root
    ns = namespace(root)
    reactions = findall("//x:reaction", root, ["x"=>ns])
    for reaction in reactions
        kineticlaws = [node for node in eachelement(reaction) if nodename(node) == "kineticLaw"]
        if length(kineticlaws) > 1
            @error("SBML reactions with more than one kinetic law are not supported.")
        end
        kineticlaw = kineticlaws[1]
        locallistsofparameters = [node for node in eachelement(kineticlaw) if occursin("Parameters",nodename(node))]
        locallistofparameters = [node for list in locallistsofparameters for node in eachelement(list)]
        for par in locallistofparameters
            if !startswith(par["id"], reaction["id"]*"_")
                _substitute_par!(kineticlaw, par["id"], reaction["id"])
                par["id"] = reaction["id"]*"_"*par["id"]
                par["name"] = reaction["id"]*"_"*par["name"]
            end
        end
    end
    doc
end

function _substitute_par!(node, par, prefix)
    if typeof(node) != Nothing
        for element in eachelement(node)
            _substitute_par!(element, par, prefix)
        end
        if strip(node.name) == "ci" && strip(node.content) == par  # && !startswith(strip(node.content), prefix*"_")
            node.content = prefix*"_"*par
        end
    end
end

""" ReactionSystem constructor """
function ModelingToolkit.ReactionSystem(sbmlmodel::SbmlModel)
    rxs = [Reaction(reac...; only_use_rate=true) for reac in sbmlmodel.reactions]
    t = Variable(:t)
    species = [spec.first for spec in sbmlmodel.species]
    pc = vcat(sbmlmodel.parameters,sbmlmodel.compartments)
    params = [par.first for par in pc]
    ReactionSystem(rxs,t,species,params)
end

""" ReactionSystem constructor """
function ModelingToolkit.ReactionSystem(sbmldocument::EzXML.Document)
    sbmlmodel = SbmlModel(sbmldocument)
    ReactionSystem(sbmlmodel)
end

""" ReactionSystem constructor """
function ModelingToolkit.ReactionSystem(sbmlfile::String)
    sbmlmodel = SbmlModel(sbmlfile)
    ReactionSystem(sbmlmodel)
end

""" ODESystem constructor """
function ModelingToolkit.ODESystem(sbmlmodel::SbmlModel)
    rs = ReactionSystem(sbmlmodel)
    u0map = [s.first => s.second[1] for s in sbmlmodel.species]
    parammap = vcat(sbmlmodel.parameters,sbmlmodel.compartments)
    defaults = vcat(u0map, parammap)
    println("System constructor")
    println(u0map)
    println(parammap)
    println(defaults)
    convert(ODESystem, rs, defaults=defaults)
end

""" ODESystem constructor """  # Todo: Why does this feed in a different SBML mdoel
function ModelingToolkit.ODESystem(sbmldocument::EzXML.Document)
    sbmlmodel = SbmlModel(sbmldocument)
    ODESystem(sbmlmodel)
end

""" ODESystem constructor """
function ModelingToolkit.ODESystem(sbmlfile::String)
    sbmlmodel = SbmlModel(sbmlfile)
    ODESystem(sbmlmodel)
end

""" ODEProblem constructor """
function ModelingToolkit.ODEProblem(sbmlmodel::SbmlModel,tspan)  # Todo: add u0 and parameters argument
    odesys = ODESystem(sbmlmodel)
    ODEProblem(odesys, [], tspan)
end

""" ODEProblem constructor """
function ModelingToolkit.ODEProblem(sbmldocument::EzXML.Document,tspan)  # Todo: add u0 and parameters argument
    odesys = ODESystem(sbmldocument)
    ODEProblem(odesys, [], tspan)
end

""" ODEProblem constructor """
function ModelingToolkit.ODEProblem(sbmlfile::String,tspan)  # Todo: add u0 and parameters argument
    odesys = ODESystem(sbmlfile)
    ODEProblem(odesys, [], tspan)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ a9afdc92-61b9-11eb-1ac4-d1b9cba169ff
using DrWatson

# ╔═╡ 98cb1be8-61bc-11eb-3652-35d883d64e97
begin
	using Pkg
end

# ╔═╡ 0d5893bc-61ba-11eb-21d6-2167c33bdeec
begin
	using ModelingToolkit, Unitful, DifferentialEquations, Optim
	using Latexify
	using BenchmarkTools, LinearAlgebra, ForwardDiff, Zygote
	using Quadrature, StaticArrays
	using Test
	using SymbolicUtils
	using SymbolicUtils: to_symbolic, substitute
	using Plots
end

# ╔═╡ 4913b034-633c-11eb-2109-5f2a85b3632f
using SICMModelingToolkit

# ╔═╡ f65c0428-61b9-11eb-0262-6dee6a26c858
DrWatson.@quickactivate "SICMModelingToolkit"

# ╔═╡ f00943f2-61db-11eb-025f-27e4a44b10f8
MTderivative(O, v) = ModelingToolkit.derivative(O, v; simplify=true)

# ╔═╡ 50a2ab6a-61bc-11eb-3933-890b4a8dfb64
rdot(v1, v2) = sum(v1 .* v2)

# ╔═╡ ef680ac4-6338-11eb-0ae4-83ed341a561e
@macroexpand @variables u(t)

# ╔═╡ 9c77aa44-6339-11eb-17b2-4fc2734cdcc6
@variables u(t)

# ╔═╡ 3b72e1d2-6339-11eb-3242-65b0429b5198
@macroexpand @parameters t2

# ╔═╡ 456b0be4-6339-11eb-39d8-17e640984fd1
t3 = Num(Sym(:t3))

# ╔═╡ 54322932-633c-11eb-0b94-b3a11c262da9
SICMModelingToolkit.generate_generic_localtuple(3)

# ╔═╡ 39738640-633c-11eb-17d0-c9dac4a1f5bb
"q1"*"_t"

# ╔═╡ f104c252-633b-11eb-233e-6d728a946be0
qlength = 3

# ╔═╡ f690e782-633b-11eb-1f36-71ab4e952e40
qnames = [Symbol("q$(i)") for i in 1:qlength]

# ╔═╡ 5c5d289e-6339-11eb-318f-15341be56c4b
u1 = (Num)(((Sym){(SymbolicUtils.FnType){NTuple{1, Any}, Real}}(:u1))((ModelingToolkit.value)(t)))

# ╔═╡ 87cdd570-6339-11eb-2f0d-e5e1ac662fab
typeof(u1)

# ╔═╡ 949a44a8-6339-11eb-3e2c-ddfbed92f6af
dump(u)

# ╔═╡ b004abe8-6339-11eb-3f8b-aba28292a88e
dump(u1)

# ╔═╡ b4d037e6-6339-11eb-12b4-25efa11e3aba
println("clear")

# ╔═╡ 4de82ab6-6339-11eb-16d4-b51a0aec0802
dump(t3)

# ╔═╡ 79a3f678-61cf-11eb-1d7a-31224592dc85
###### 1.4

# ╔═╡ 1a88e708-61ba-11eb-394a-fff715bf6808
@parameters t m k

# ╔═╡ 0ffa5e60-61bb-11eb-05e1-0d2a6f66ed62
@variables x(t) y(t) z(t)

# ╔═╡ 78d740de-61dd-11eb-1863-cb24f75a575a
@variables x_t(t) y_t(t) z_t(t)

# ╔═╡ 145d5854-61bb-11eb-18db-f74e83fdde06
D = Differential(t)

# ╔═╡ 1eb34a04-61bb-11eb-0c4b-8d3e4a85d3e9
q = @SVector [x, y, z]

# ╔═╡ 96d39b7a-61bb-11eb-2170-1178df4203f0
qv = @SVector [x_t, y_t, z_t]

# ╔═╡ 4936584c-61ba-11eb-1d82-11d5fd37d2e1
L_free_particle(t,x,v) = 1//2 .* m .* rdot(v,v)

# ╔═╡ be3db274-61bb-11eb-3dc1-3bb68cc3caf8
latexify(L_free_particle(t,q,qv))

# ╔═╡ c2f96f38-61bb-11eb-1d0b-bd366c75515e
q .+ q

# ╔═╡ 69333044-61bc-11eb-0784-7d2954fff119
qeqs = q .~ @SVector [4t + 7, 3t + 5, 2t + 1]

# ╔═╡ 37011c9e-61be-11eb-0562-3b7a61c8cd5b
qns = NonlinearSystem(qeqs, q, @SVector [t])

# ╔═╡ bc4ae90c-61be-11eb-1634-45317900a809
qfunc = eval(generate_function(qns)[2])

# ╔═╡ 48367b52-61bf-11eb-3cef-e3067a350e8b
begin
	du = zeros(3)
	qfunc(du, q, [-1.0])
end

# ╔═╡ c199edee-61bf-11eb-3e76-d7c1c181f14b
du

# ╔═╡ 2c2bb89c-61c0-11eb-1526-9fdc91a344d9
ModelingToolkit.substitute(qeqs, Dict(t=>1, x=>5t + 7))

# ╔═╡ d12126f4-61c1-11eb-3a8b-215ff2d0aae4
map(qeq->substitute(qeq, Dict(x=>5t+7, y=>3t+5, z=>2t+1)), q)

# ╔═╡ a8b2d568-61c2-11eb-300f-757381afee6a
test_path = @SVector [4t + 7, 3t + 5, 2t + 1]

# ╔═╡ 61393e6e-61dd-11eb-1d03-950972079936
@variables x_t(t)

# ╔═╡ 31f036b6-61c2-11eb-3afa-35bdfaffbad4
function Γ(qvars, coordfuncs; numD=2)
	subdict = Dict(qvars .=> coordfuncs)
	mapsubs(vals) = map(qvar->expand_derivatives(substitute(qvar, subdict)), vals)
	subq = mapsubs(qvars)
	if numD > 0
		derivativeqs = [mapsubs((D^i).(qvars)) for i in 1:numD]
		[[t]; [subq]; derivativeqs]
	else
		[[t]; [subq]]
	end
end

# ╔═╡ 0467956e-61dd-11eb-2604-81571591b6c3
function Γ_vars(qvars; numD=2)
	subdict = Dict(qvars .=> coordfuncs)
	mapsubs(vals) = map(qvar->expand_derivatives(substitute(qvar, subdict)), vals)
	subq = mapsubs(qvars)
	if numD > 0
		derivativeqs = [mapsubs((D^i).(qvars)) for i in 1:numD]
		[[t]; [subq]; derivativeqs]
	else
		[[t]; [subq]]
	end
end

# ╔═╡ 7c4f9a38-61c2-11eb-0d94-bf3bf88b0875
Γ(q, test_path; numD=2)[4]

# ╔═╡ e2108a08-61c2-11eb-0247-d52e97b2ec87
function compose_L_Γ(L_i)
	function composed(qvars, coordfuncs; numD=2)
		Γ_val = Γ(qvars, coordfuncs; numD=numD)
		L_i(Γ_val...)
	end
end

# ╔═╡ eaf1913c-61c5-11eb-04ee-29014e6d0f8b
L_free_composed_test_path = compose_L_Γ(L_free_particle)(q, test_path; numD=1)

# ╔═╡ 3c424016-61c8-11eb-1c0a-d3112cf4ff55
L_free_eval(1.0, 3.0)

# ╔═╡ 947249ce-61c9-11eb-23f3-3d246bd8dc0b
function Lagrangian_action(L_i, path_func, t_1, t_2, qvars, pvars, pvals; numD=2)
	L_composed = compose_L_Γ(L_i)(qvars, path_func; numD=numD)
	L_func = build_function(L_composed, [t], pvars; expression=Val{false})
	quadprob = QuadratureProblem(L_func, t_1, t_2, pvals)
	solve(quadprob, QuadGKJL(); abstol=1e-7).u
end

# ╔═╡ 0c470b92-61ca-11eb-3cde-d5112faac779
func = Lagrangian_action(L_free_particle, test_path, 0.0, 10.0, q, [m], 3.0; numD=1)

# ╔═╡ 6737dc14-61cc-11eb-3ee1-85ccdd31619c
function make_η(ν, t_1, t_2)
	(t .- t_1) .* (t .- t_2) .* ν
end

# ╔═╡ 8b97dbfe-61cc-11eb-1545-9bc13cdd6614
example_ν = @SVector [sin(t), cos(t), t^2]

# ╔═╡ a4174b44-61cc-11eb-06cf-873f32b00f21
make_η(example_ν, 0.0, 10.0)

# ╔═╡ 6a5de2d6-61cb-11eb-0d7b-538845bef060
function varied_free_particle_action(mass, path, ν, t_1, t_2)
	η = make_η(ν, t_1, t_2)
	function per_eps(ϵ)
		perturbed_path = path .+ ϵ .* η
		Lagrangian_action(L_free_particle, perturbed_path, t_1, t_2, q, [m], mass; numD=1)
	end
end

# ╔═╡ 03cefa1a-61cd-11eb-2cba-2f7b8a9f4fc0
varied_free_action_per_ϵ = varied_free_particle_action(3.0, test_path, example_ν, 0.0, 10.0)

# ╔═╡ 37ec4df4-61cd-11eb-295f-e904a0abc895
varied_free_action_per_ϵ(0.001)

# ╔═╡ 47cbd212-61cd-11eb-1899-274b2e205ead
Optim.optimize(varied_free_action_per_ϵ, -2.0, 1.0)

# ╔═╡ ef791fe8-61ce-11eb-156e-b5b67a0e7887
#### 1.5

# ╔═╡ 39ae9600-61d0-11eb-0aed-63749ad212d0
Differential(t)(x)

# ╔═╡ 76e0551c-61cf-11eb-26a1-7522376db177
expand_derivatives(Differential(D(x))(L_free_particle(Γ(q, q; numD=1)...)))

# ╔═╡ 079bbdd8-61d2-11eb-0092-fb7654b3b0ae
Differential(x)(L_free_particle(Γ(q, q; numD=1)...))

# ╔═╡ e66f3d66-61d0-11eb-3dbb-03372172ba86
L_free_particle(Γ(q, q; numD=1)...)

# ╔═╡ eb0fb58a-61df-11eb-1f5b-c731c5ad7ece
function find_multiplies_to_derivative(expression)
	#@show expression
	#@show typeof(expression.val)
	#@show ModelingToolkit.isdifferential(expression.val)
	diffs = []
	premults = []
	if typeof(expression) <: SymbolicUtils.Mul
		#@show "mul!!"
		coeff = expression.coeff
		push!(premults, coeff)
		for key in keys(expression.dict)
			if ModelingToolkit.isdifferential(key)
				#@show "isa diff"
				push!(diffs, key)
			else
				#@show "isnota diff"
				push!(premults, key)
			end
			#@show key
		end
		#recursed_mul_premult, recursed_mul_deriv = recursively_find_multiplies_to_derivative(
	elseif ModelingToolkit.isdifferential(expression)
		push!(diffs, expression)
	end
	diffs, premults
end

# ╔═╡ 13a053bc-61e0-11eb-0802-9d678d04e727
diffs, premults = find_multiplies_to_derivative((m * D(x)).val)

# ╔═╡ b3be1294-61e3-11eb-298a-5f6e59e7b0df
length(premults)

# ╔═╡ 172bc650-61e4-11eb-3514-156630da8818
dump(prod(premults))

# ╔═╡ 1232c1e8-61e5-11eb-1a93-e74e64faab02
typeof(prod(premults))

# ╔═╡ b0066b96-61e4-11eb-1791-0d5c9a8032e0
typeof(Num(diffs[1]))

# ╔═╡ 2bc7d170-61e5-11eb-22e7-198087098446
typeof(D(x))

# ╔═╡ c54e959a-61e0-11eb-12ac-e35a69d26a18
keys((m * x).val.dict)

# ╔═╡ 15ca6e50-61d1-11eb-1f59-1b814635c335
function Lagrange_equations(L_i, qvars, qdotvars)
	symbolic_L = L_i(t, qvars, qdotvars)
	equations = []
	for (qvar, qdotvar) in zip(qvars, qdotvars)
		lhs_premult = MTderivative(MTderivative(symbolic_L, qdotvar), t)
		lhs_diffs, lhs_mults = find_multiplies_to_derivative(lhs_premult.val)
		rhs = MTderivative(symbolic_L, qvar)
		lhs = lhs_premult
		if length(lhs_mults) > 0
			lhs_multiplied = prod(Num.(lhs_mults))
			rhs = rhs / lhs_multiplied
			lhs = Num(lhs_diffs[1])
		end
		push!(equations, lhs ~ rhs)
		push!(equations, D(qvar) ~ qdotvar)
	end
	equations
end

# ╔═╡ a4ae8bd0-62ab-11eb-1f76-2596a3d649bb

# ╔═╡ 61d6f634-61da-11eb-3621-4bae15be5cfe
length(x)

# ╔═╡ 45a8e8d6-61d1-11eb-2e54-97bd0832cf66
lagfree = Lagrange_equations(L_free_particle, q, qv)

# ╔═╡ c6c42192-61d1-11eb-0df7-a1300ccb1fa0
lagfree[3].lhs

# ╔═╡ 2a52fcf6-61d2-11eb-3632-217ef932b5c1
lagrangian_eqs_sys = ODESystem(lagfree, t, q, qv, [m])

# ╔═╡ 63350bcc-61d2-11eb-27c1-5b88ddbb29de
Lagrangian_eqs_sys.eq

# ╔═╡ b61ca1ce-61d2-11eb-2bb9-8be89a023b1f
lagrangian_eqs_sys.eqs

# ╔═╡ 1d6d2c86-61d3-11eb-1752-4df7966d2898
@parameters k

# ╔═╡ 10953204-61d3-11eb-23ee-bd299a2e470a
L_harmonic_spring(tvar, qvar, qdotvar) = 1//2 * m * qdotvar^2 - 1//2 * k * qvar^2

# ╔═╡ 3e7fe27c-61d3-11eb-160e-39eeba93ac65
L_harmonic_spring(t, x, D(x))

# ╔═╡ 551f58ac-61d3-11eb-0fdf-79af470d57bd
lagharm = Lagrange_equations(L_harmonic_spring, x, x_t)

# ╔═╡ 8496a978-61d3-11eb-3ef9-35e8f3f18f6a
lagharm

# ╔═╡ 9d53a218-61d3-11eb-3395-bf43539f832f
harmonic_sys = ODESystem(lagharm, t)

# ╔═╡ e56035c2-61dc-11eb-20f0-b1e5e6664dd4
harmonic_sys.eqs

# ╔═╡ 90da45a6-61e1-11eb-2f00-11bba22507d7
harmonic_sys.ps[1].name

# ╔═╡ 26dd26f4-61e2-11eb-3f53-15c38b38e1c8
harmonic_sys.states[2].f.name

# ╔═╡ 4d933d2c-61e2-11eb-0433-5bcd8219b181

# ╔═╡ d7a1a312-61d4-11eb-332a-7124267dd164
begin
	tspan = (0.0, π)
	u0 = [1.0, 0.0]
	p0 = [1.0, 3.0]
end

# ╔═╡ be6c49f4-61d4-11eb-018d-079c83d1ad64
harmonic_prob = ODEProblem(harmonic_sys, [x=>u0[1], x_t=>u0[2]], tspan, [k=>p0[1], m=>p0[2]])

# ╔═╡ d0d1dab4-61d9-11eb-3a8e-7b5d67dfaf22
harmonic_sol = solve(harmonic_prob, Tsit5())

# ╔═╡ e539f43c-61d9-11eb-0453-9d96c45fcc95
plot(harmonic_sol; dpi=350)

# ╔═╡ ac536e02-61d7-11eb-1521-193caca80dbe
testeqs = [D(x) ~ 1]

# ╔═╡ bb348136-61d7-11eb-36d8-d7b6a9f28c26
testsys = ODESystem(testeqs, t, [x], [])

# ╔═╡ cabc1f6a-61d7-11eb-3bc1-a3d75e1f8058
testfunc = ODEFunction(testsys, [x], [])

# ╔═╡ f4b54fe4-61d7-11eb-19f1-3d9a86d5e46d
testprob = ODEProblem(testsys, [0.0], (0.0, 2.0))

# ╔═╡ 43cbdb98-61d8-11eb-2c8b-fded5290943d
testsol = solve(testprob, Tsit5())

# ╔═╡ 4a27a706-61d8-11eb-3d04-c5d7323a7249
plot(testsol; dpi=350)

# ╔═╡ 82cf65a0-61d8-11eb-24b3-1b1beb77e372
testsecondordereqs = [D(D(x)) ~ 1]

# ╔═╡ 9d33eee6-61d8-11eb-2a45-41a7849fc0ca
testsecondordersys = ODESystem(testsecondordereqs)

# ╔═╡ b49bc720-61d8-11eb-2ff5-099ee4ac56c0
testfosys = ode_order_lowering(testsecondordersys)

# ╔═╡ 75b20a8e-61d9-11eb-1bbf-6331e9accaf1
testfosys.eqs

# ╔═╡ cb50bb58-61d8-11eb-210d-0bbadabc5508
testfosys.states

# ╔═╡ dcac2160-61d8-11eb-1d82-b934f5d2eae7
testfoprob = ODEProblem(testfosys, [0.0, 0.0], (0.0, 4.0))

# ╔═╡ ef86b534-61d8-11eb-1c68-630277726bf4
testfosol = solve(testfoprob, Tsit5())

# ╔═╡ f9cadcc6-61d8-11eb-3fd1-374577903088
plot(testfosol; dpi=350)

# ╔═╡ 41f7fcac-61dd-11eb-37b5-dd1b334525e5
### kepler problem
### k = μ, m = m, x = η, y = ζ

# ╔═╡ 9f9ae924-61e6-11eb-2150-2ffd405f3a1e
L_kepler(t, q, qv) = 1//2 * m * rdot(qv, qv) + k / sqrt(rdot(q, q))

# ╔═╡ f7a3b9c2-61e6-11eb-2bba-f7d1d6fe1658
L_kepler(t, [x, y], [x_t, y_t])

# ╔═╡ 045eda50-61e7-11eb-136b-47bb20bf6c6e
lagkepler = Lagrange_equations(L_kepler, [x, y], [x_t, y_t])

# ╔═╡ 2b3bcca8-61e9-11eb-204f-4d0e844680cd
kepler_sys = ODESystem(lagkepler, t)

# ╔═╡ 76155afa-6211-11eb-2d38-0f88ece4aabf
kepler_ode_func = ODEFunction(kepler_sys)

# ╔═╡ de9457d4-6211-11eb-020c-5d19dcaf8106
begin
	u0_kep = [-1.0, 0.0, 0.0, 1.0]
	du_kep = similar(u0_kep)
	tspan_kep = (0.0, 5.0)
	p_kep = [10.0, 3.0]
end

# ╔═╡ 229ec7d8-6213-11eb-261d-4739e7562178
kepler_ode_func(du_kep, u0_kep, p_kep, 0.0)

# ╔═╡ 3cba55e2-6213-11eb-2570-cf445b02d7c5
du_kep

# ╔═╡ 22c98e36-61e7-11eb-0389-2fcfe5f3bd11
kepler_prob = ODEProblem(kepler_ode_func, u0_kep, tspan_kep, p_kep)

# ╔═╡ 3c58cd92-61e9-11eb-34ec-c1b507fdffc2
kepler_sol = solve(kepler_prob, Vern8())

# ╔═╡ 434fd42e-61e9-11eb-3178-e76dc841018f
kepplot = plot(kepler_sol; dpi=350)

# ╔═╡ e291cf0c-6210-11eb-2fca-6b83b66c9de3
kepphaseplot = plot(kepler_sol; vars=(2, 4), dpi=350)

# ╔═╡ 11dffbf8-6211-11eb-382b-dba7dfe9ed2a
kepphasetimeplot = plot(kepler_sol; vars=(0, 2, 4), dpi=350)

# ╔═╡ b1430fd4-620f-11eb-1910-15102c12f52f
save("kepplot.png", kepplot)

# ╔═╡ 18127514-6211-11eb-07f8-33293bb5683b
save("kepphaseplot.png", kepphaseplot)

# ╔═╡ 24b3de8e-6211-11eb-2901-e30917142b8d
save("kepphasetimeplot.png", kepphasetimeplot)

# ╔═╡ 3823af7e-61fa-11eb-2095-1f864cf71459
#### 1.6.1 Coordinate Transformations

# ╔═╡ 4402830e-61fa-11eb-2c55-354380e16398
function F_to_C(F)
	function F_to_C_inner(locals)
		(t, F(locals), MTderivative(F, locals[2])(locals) .+ MTderivative(F, locals[3])(locals) .* locals[3])
	end
end

# ╔═╡ 0455cb86-61fb-11eb-2ccf-7786a457e9e0
function polar_to_rectangular(locals)
	r, ϕ = locals[2]
	x = r * cos(ϕ)
	y = r * sin(ϕ)

	(x, y)
end

# ╔═╡ ed46d722-61fb-11eb-1995-833ff65293be
@variables r(t) ϕ(t)

# ╔═╡ 0aeed2de-61fc-11eb-2955-b5a0b54abe00
@variables r_t(t) ϕ_t(t)

# ╔═╡ f41cb288-61fb-11eb-07a4-31db6009cfc7
locals = (t, (r, ϕ), (r_t, ϕ_t))

# ╔═╡ 1841a236-61fc-11eb-1753-d392c6a20d57
polar_to_rectangular(locals)

# ╔═╡ 539f0c16-61fe-11eb-2141-fbae47ab042e
prl1 = polar_to_rectangular(locals)[1]

# ╔═╡ c171da40-61fe-11eb-161f-e3d6cbea4ab2
prl2 = polar_to_rectangular(locals)[2]

# ╔═╡ 5bcb895c-61fe-11eb-31bf-e51353efd58d
simplify(MTderivative(prl1, t) + MTderivative(prl1, r) * D(r) + MTderivative(prl1, ϕ) * D(ϕ))

# ╔═╡ c036a7a0-61fe-11eb-0c48-97d99436034b
simplify(MTderivative(prl2, t) + MTderivative(prl2, r) * D(r) + MTderivative(prl2, ϕ) * D(ϕ))

# ╔═╡ 52aa43a6-61fc-11eb-1501-4ff921328b78
expanded_x = simplify(MTderivative(polar_to_rectangular(locals)[2], locals[1]) + MTderivative(polar_to_rectangular(locals)[2], locals[2][1]) * locals[3][1])

# ╔═╡ 067c243a-61fd-11eb-131f-c96e8cb2aae5
derivativerrule = @acrule (D(~r)) => ~r_t

# ╔═╡ 8701b8ea-61fd-11eb-3b03-35ba2013864d
derivativerrule(expanded_x)

# ╔═╡ 77f6c2ec-61fc-11eb-1c5b-436ba55225be
polar_to_rectangular(locals)[1]

# ╔═╡ 28bdcd58-61fc-11eb-22fc-99210bba2842
F_to_C(polar_to_rectangular)(locals)

# ╔═╡ Cell order:
# ╠═a9afdc92-61b9-11eb-1ac4-d1b9cba169ff
# ╠═f65c0428-61b9-11eb-0262-6dee6a26c858
# ╠═98cb1be8-61bc-11eb-3652-35d883d64e97
# ╠═0d5893bc-61ba-11eb-21d6-2167c33bdeec
# ╠═4913b034-633c-11eb-2109-5f2a85b3632f
# ╠═f00943f2-61db-11eb-025f-27e4a44b10f8
# ╠═50a2ab6a-61bc-11eb-3933-890b4a8dfb64
# ╠═ef680ac4-6338-11eb-0ae4-83ed341a561e
# ╠═9c77aa44-6339-11eb-17b2-4fc2734cdcc6
# ╠═3b72e1d2-6339-11eb-3242-65b0429b5198
# ╠═456b0be4-6339-11eb-39d8-17e640984fd1
# ╠═54322932-633c-11eb-0b94-b3a11c262da9
# ╠═39738640-633c-11eb-17d0-c9dac4a1f5bb
# ╠═f104c252-633b-11eb-233e-6d728a946be0
# ╠═f690e782-633b-11eb-1f36-71ab4e952e40
# ╠═5c5d289e-6339-11eb-318f-15341be56c4b
# ╠═87cdd570-6339-11eb-2f0d-e5e1ac662fab
# ╠═949a44a8-6339-11eb-3e2c-ddfbed92f6af
# ╠═b004abe8-6339-11eb-3f8b-aba28292a88e
# ╠═b4d037e6-6339-11eb-12b4-25efa11e3aba
# ╠═4de82ab6-6339-11eb-16d4-b51a0aec0802
# ╠═79a3f678-61cf-11eb-1d7a-31224592dc85
# ╠═1a88e708-61ba-11eb-394a-fff715bf6808
# ╠═0ffa5e60-61bb-11eb-05e1-0d2a6f66ed62
# ╠═78d740de-61dd-11eb-1863-cb24f75a575a
# ╠═145d5854-61bb-11eb-18db-f74e83fdde06
# ╠═1eb34a04-61bb-11eb-0c4b-8d3e4a85d3e9
# ╠═96d39b7a-61bb-11eb-2170-1178df4203f0
# ╠═4936584c-61ba-11eb-1d82-11d5fd37d2e1
# ╠═be3db274-61bb-11eb-3dc1-3bb68cc3caf8
# ╠═c2f96f38-61bb-11eb-1d0b-bd366c75515e
# ╠═69333044-61bc-11eb-0784-7d2954fff119
# ╠═37011c9e-61be-11eb-0562-3b7a61c8cd5b
# ╠═bc4ae90c-61be-11eb-1634-45317900a809
# ╠═48367b52-61bf-11eb-3cef-e3067a350e8b
# ╠═c199edee-61bf-11eb-3e76-d7c1c181f14b
# ╠═2c2bb89c-61c0-11eb-1526-9fdc91a344d9
# ╠═d12126f4-61c1-11eb-3a8b-215ff2d0aae4
# ╠═a8b2d568-61c2-11eb-300f-757381afee6a
# ╠═61393e6e-61dd-11eb-1d03-950972079936
# ╠═31f036b6-61c2-11eb-3afa-35bdfaffbad4
# ╠═0467956e-61dd-11eb-2604-81571591b6c3
# ╠═7c4f9a38-61c2-11eb-0d94-bf3bf88b0875
# ╠═e2108a08-61c2-11eb-0247-d52e97b2ec87
# ╠═eaf1913c-61c5-11eb-04ee-29014e6d0f8b
# ╠═3c424016-61c8-11eb-1c0a-d3112cf4ff55
# ╠═947249ce-61c9-11eb-23f3-3d246bd8dc0b
# ╠═0c470b92-61ca-11eb-3cde-d5112faac779
# ╠═6737dc14-61cc-11eb-3ee1-85ccdd31619c
# ╠═8b97dbfe-61cc-11eb-1545-9bc13cdd6614
# ╠═a4174b44-61cc-11eb-06cf-873f32b00f21
# ╠═6a5de2d6-61cb-11eb-0d7b-538845bef060
# ╠═03cefa1a-61cd-11eb-2cba-2f7b8a9f4fc0
# ╠═37ec4df4-61cd-11eb-295f-e904a0abc895
# ╠═47cbd212-61cd-11eb-1899-274b2e205ead
# ╠═ef791fe8-61ce-11eb-156e-b5b67a0e7887
# ╠═39ae9600-61d0-11eb-0aed-63749ad212d0
# ╠═76e0551c-61cf-11eb-26a1-7522376db177
# ╠═079bbdd8-61d2-11eb-0092-fb7654b3b0ae
# ╠═e66f3d66-61d0-11eb-3dbb-03372172ba86
# ╠═eb0fb58a-61df-11eb-1f5b-c731c5ad7ece
# ╠═13a053bc-61e0-11eb-0802-9d678d04e727
# ╠═b3be1294-61e3-11eb-298a-5f6e59e7b0df
# ╠═172bc650-61e4-11eb-3514-156630da8818
# ╠═1232c1e8-61e5-11eb-1a93-e74e64faab02
# ╠═b0066b96-61e4-11eb-1791-0d5c9a8032e0
# ╠═2bc7d170-61e5-11eb-22e7-198087098446
# ╠═c54e959a-61e0-11eb-12ac-e35a69d26a18
# ╠═15ca6e50-61d1-11eb-1f59-1b814635c335
# ╠═a4ae8bd0-62ab-11eb-1f76-2596a3d649bb
# ╠═61d6f634-61da-11eb-3621-4bae15be5cfe
# ╠═45a8e8d6-61d1-11eb-2e54-97bd0832cf66
# ╠═c6c42192-61d1-11eb-0df7-a1300ccb1fa0
# ╠═2a52fcf6-61d2-11eb-3632-217ef932b5c1
# ╠═63350bcc-61d2-11eb-27c1-5b88ddbb29de
# ╠═b61ca1ce-61d2-11eb-2bb9-8be89a023b1f
# ╠═1d6d2c86-61d3-11eb-1752-4df7966d2898
# ╠═10953204-61d3-11eb-23ee-bd299a2e470a
# ╠═3e7fe27c-61d3-11eb-160e-39eeba93ac65
# ╠═551f58ac-61d3-11eb-0fdf-79af470d57bd
# ╠═8496a978-61d3-11eb-3ef9-35e8f3f18f6a
# ╠═9d53a218-61d3-11eb-3395-bf43539f832f
# ╠═e56035c2-61dc-11eb-20f0-b1e5e6664dd4
# ╠═90da45a6-61e1-11eb-2f00-11bba22507d7
# ╠═26dd26f4-61e2-11eb-3f53-15c38b38e1c8
# ╠═4d933d2c-61e2-11eb-0433-5bcd8219b181
# ╠═d7a1a312-61d4-11eb-332a-7124267dd164
# ╠═be6c49f4-61d4-11eb-018d-079c83d1ad64
# ╠═d0d1dab4-61d9-11eb-3a8e-7b5d67dfaf22
# ╠═e539f43c-61d9-11eb-0453-9d96c45fcc95
# ╠═ac536e02-61d7-11eb-1521-193caca80dbe
# ╠═bb348136-61d7-11eb-36d8-d7b6a9f28c26
# ╠═cabc1f6a-61d7-11eb-3bc1-a3d75e1f8058
# ╠═f4b54fe4-61d7-11eb-19f1-3d9a86d5e46d
# ╠═43cbdb98-61d8-11eb-2c8b-fded5290943d
# ╠═4a27a706-61d8-11eb-3d04-c5d7323a7249
# ╠═82cf65a0-61d8-11eb-24b3-1b1beb77e372
# ╠═9d33eee6-61d8-11eb-2a45-41a7849fc0ca
# ╠═b49bc720-61d8-11eb-2ff5-099ee4ac56c0
# ╠═75b20a8e-61d9-11eb-1bbf-6331e9accaf1
# ╠═cb50bb58-61d8-11eb-210d-0bbadabc5508
# ╠═dcac2160-61d8-11eb-1d82-b934f5d2eae7
# ╠═ef86b534-61d8-11eb-1c68-630277726bf4
# ╠═f9cadcc6-61d8-11eb-3fd1-374577903088
# ╠═41f7fcac-61dd-11eb-37b5-dd1b334525e5
# ╠═9f9ae924-61e6-11eb-2150-2ffd405f3a1e
# ╠═f7a3b9c2-61e6-11eb-2bba-f7d1d6fe1658
# ╠═045eda50-61e7-11eb-136b-47bb20bf6c6e
# ╠═2b3bcca8-61e9-11eb-204f-4d0e844680cd
# ╠═76155afa-6211-11eb-2d38-0f88ece4aabf
# ╠═229ec7d8-6213-11eb-261d-4739e7562178
# ╠═3cba55e2-6213-11eb-2570-cf445b02d7c5
# ╠═de9457d4-6211-11eb-020c-5d19dcaf8106
# ╠═22c98e36-61e7-11eb-0389-2fcfe5f3bd11
# ╠═3c58cd92-61e9-11eb-34ec-c1b507fdffc2
# ╠═434fd42e-61e9-11eb-3178-e76dc841018f
# ╠═e291cf0c-6210-11eb-2fca-6b83b66c9de3
# ╠═11dffbf8-6211-11eb-382b-dba7dfe9ed2a
# ╠═b1430fd4-620f-11eb-1910-15102c12f52f
# ╠═18127514-6211-11eb-07f8-33293bb5683b
# ╠═24b3de8e-6211-11eb-2901-e30917142b8d
# ╠═3823af7e-61fa-11eb-2095-1f864cf71459
# ╠═4402830e-61fa-11eb-2c55-354380e16398
# ╠═0455cb86-61fb-11eb-2ccf-7786a457e9e0
# ╠═ed46d722-61fb-11eb-1995-833ff65293be
# ╠═0aeed2de-61fc-11eb-2955-b5a0b54abe00
# ╠═f41cb288-61fb-11eb-07a4-31db6009cfc7
# ╠═1841a236-61fc-11eb-1753-d392c6a20d57
# ╠═539f0c16-61fe-11eb-2141-fbae47ab042e
# ╠═c171da40-61fe-11eb-161f-e3d6cbea4ab2
# ╠═5bcb895c-61fe-11eb-31bf-e51353efd58d
# ╠═c036a7a0-61fe-11eb-0c48-97d99436034b
# ╠═52aa43a6-61fc-11eb-1501-4ff921328b78
# ╠═067c243a-61fd-11eb-131f-c96e8cb2aae5
# ╠═8701b8ea-61fd-11eb-3b03-35ba2013864d
# ╠═77f6c2ec-61fc-11eb-1c5b-436ba55225be
# ╠═28bdcd58-61fc-11eb-22fc-99210bba2842

max_stars_repo_name: zoemcc/SICMModelingToolkit

max_issues_repo_name: zoemcc/SICMModelingToolkit

max_forks_repo_name: zoemcc/SICMModelingToolkit


text:
function set_linearized_b!(b, q′, d::Domain)
    c = CartesianIndices(d)
    l = LinearIndices(d)
    for I = c
        b[l[I]] = q′[I]
    end
    ol = last(l)
    for I = c
        b[ol+l[I]] = 0
    end
end

function set_bψ!(bψ, ψ, ϕ, q, d::Domain)
    c = CartesianIndices(d)
    l = LinearIndices(d)
    for I = c
        bψ[l[I]] = (
            q[I] + ∂x²(I,ϕ,d) + ∂y²(I,ϕ,d) - ∂z²(I,ϕ,d)
          - 2*(∂x²(I,ψ,d)*∂y²(I,ψ,d) - ∂xy(I,ψ,d)^2)
          + ∂xz(I,ψ,d)*∂xz(I,ϕ,d) + ∂yz(I,ψ,d)*∂yz(I,ϕ,d)
        )
    end
    return bψ
end

function set_bϕ!(bϕ, ψ, ϕ, q, d::Domain)
    c = CartesianIndices(d)
    l = LinearIndices(d)
    for I = c
        bϕ[l[I]] = (
            q[I] + ∂x²(I,ψ,d) + ∂y²(I,ψ,d)
          + 2*(∂x²(I,ψ,d)*∂y²(I,ψ,d) - ∂xy(I,ψ,d)^2)
          + ∂xz(I,ψ,d)*∂xz(I,ϕ,d) + ∂yz(I,ψ,d)*∂yz(I,ϕ,d)
        )
    end
    return bϕ
end

function set_cv_bψ!(bψ, q, d::Domain)
    c = CartesianIndices(d)
    l = LinearIndices(d)
    for I = c
        bψ[l[I]] = q[I] - 1
    end
    return bψ
end
function set_cv_bϕ!(bϕ, ψ, d::Domain)
    c = CartesianIndices(d)
    l = LinearIndices(d)
    for I = c
        bϕ[l[I]] = (
            ∂x²(I,ψ,d) + ∂y²(I,ψ,d)
          + 2*(∂x²(I,ψ,d)*∂y²(I,ψ,d) - (∂xy(I,ψ,d)^2))
        )
    end
    return bϕ
end


text: while (x = "41" + 1; true) end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
module F
export f

function f(arg::Vector{Float64})
    return sum(arg)
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
    AbstractGraph{V,E}

Abstract supertype for all graph types, where `V` represents node type and `E` represents edge type.
"""
abstract type AbstractGraph{V,E} end

"""
    AbstractEvolvingGraph{V,E,T} <: AbstractGraph{V,E}

Abstract supertype for all evolving graph types, where `V` represents node type, `E` represents edge type, and `T` represents timestamp type.
"""
abstract type AbstractEvolvingGraph{V,E,T} <: AbstractGraph{V,E} end

"""
    AbstractStaticGraph{V,E} <: AbstractGraph{V,E}

Abstract supertype for all static graph types, where `V` represents node type and `E` represents edge etype.
"""
abstract type AbstractStaticGraph{V,E} <: AbstractGraph{V,E} end

"""
    AbstractNode{V}

Abstract supertype for all node types, where `V` is the type of the node key.
"""
abstract type AbstractNode{V} end

"""
    node_index(v)
    node_index(g::AbstractGraph, v)

Return the index of a node `v`.
"""
node_index(v::AbstractNode) = v.index
node_index(g::AbstractGraph, v::AbstractNode) = node_index(v)

"""
    node_key(v)

Return the key of a node `v`.
"""
node_key(v::AbstractNode) = v.key

"""
    Node(index, key)
    Node{V}(key)
    Node{V}(g, key)

Constructs a Node with node index `index` of type `V` and key value `key`.
If only `key` is given, set `index` to `0`. Node(g, key) constructs a new node in graph `g`, so index is equal to `num_nodes(g) + 1`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> a = Node(1,"a")
Node(a)

julia> node_index(a)
1

julia> node_key(a)
"a"

julia> b = Node{String}("b")
Node(b)

julia> node_index(b)
0
```
"""
struct Node{V} <: AbstractNode{V}
    index::Int
    key::V
end
Node{V}(key::V) where {V} = Node(0, key)
Node{V}(g::AbstractGraph, key::V) where {V} = Node(num_nodes(g)+ 1, key)

==(v1::Node, v2::Node) = (v1.key == v2.key) && (v1.index == v2.index)
eltype{T}(::Node{T}) = T
eltype{T}(::Type{Node{T}}) = T

"""
    AttributeNode(index, key, attributes=Dict())
    AttributeNode{K, D_k, D_v}(key, attributes=Dict())
    AttributeNode{K, D_k, D_v}(g, key, attributes=Dict())

Construct an AttributeNode with index `index`, key value `key` and attributes `attributes`. If `index` is not given, set `index = 0`. `AttributeNode(g, key, attributes)` constructs a new AttributeNode for graph `g`, where `index = num_nodes(g) + 1`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> a = AttributeNode(1, "a", Dict("a" => 12))
AttributeNode(a)

julia> node_key(a)
"a"

julia> node_index(a)
1

julia> node_attributes(a)
Dict{String,Int64} with 1 entry:
  "a" => 12
```
"""
mutable struct AttributeNode{V, D_k, D_v} <: AbstractNode{V}
    index::Int
    key::V
    attributes::Dict{D_k, D_v}
end
AttributeNode{V, D_k, D_v}(key::V, attributes=Dict{D_k, D_v}()) where V where D_k where D_v = AttributeNode(0, key, attributes)
AttributeNode{V}(key::V, attributes=Dict()) where V = AttributeNode(0, key, attributes)
AttributeNode{V, D_k, D_v}(g::AbstractGraph, key::V, attributes=Dict()) where V where D_k where D_v = AttributeNode(num_nodes(g) + 1 , key, attributes)
AttributeNode{KV}(g::AbstractGraph, key::KV, attributes=Dict()) = AttributeNode{KV, Any, Any}(g, key, attributes=attributes)

"""
    node_attributes(v)
    node_attributes(v, g::AbstractGraph)

Returns the attributes of AttributeNode `v`.
"""
node_attributes(v::AttributeNode) = v.attributes
node_attributes(v::AttributeNode, g::AbstractGraph) = v.attributes

eltype{T}(::AttributeNode{T}) = T
eltype{T}(::Type{AttributeNode{T}}) = T
==(v1::AttributeNode, v2::AttributeNode) = (v1.key == v2.key &&
                                            v1.attributes == v2.attributes && v1.index == v2.index)

"""
    TimeNode(index, key, timestamp)
    TimeNode{V,T}(key, timestamp)
    TimeNode{V,T}(g, key, timestamp)

Constructs a TimeNode at timestamp `timestamp`. `TimeNode(g, key, timestamp)` constructs a TimeNode in graph `g`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> t = TimeNode(1, "t", 2018)
TimeNode(t, 2018)

julia> node_index(t)
1

julia> node_key(t)
"t"

julia> node_timestamp(t)
2018
```
"""
struct TimeNode{V,T} <: AbstractNode{V}
    index::Int
    key::V
    timestamp::T
end
TimeNode{V,T}(key, timestamp) where {V} where {T} = TimeNode(0, key, timestamp)
TimeNode{V,T}(g::AbstractGraph, key::V, timestamp::T) where {V} where {T} = TimeNode(num_active_nodes(g) + 1, key, timestamp)

"""
    node_timestamp(v)

Return the timestamp of TimeNode `v`.
"""
node_timestamp(v::TimeNode) = v.timestamp

eltype{V,T}(::TimeNode{V,T}) = (V, T)
==(v1::TimeNode, v2::TimeNode) = (v1.key == v2.key && v1.timestamp== v2.timestamp)

"""
    AbstractEdge{V}

Abstract supertype for all edge types, where `V` represents the type of the source and target nodes.
"""
abstract type AbstractEdge{V} end

"""
    source(e)

Return the source of edge `e`.
"""
source(e::AbstractEdge) = e.source

"""
    target(e)

Return the target of edge `e`.
"""
target(e::AbstractEdge) = e.target

"""
    edge_timestamp(e)

Return the timestamp of edge `e` if `e` is a time dependent edge, i.e., `TimeEdge` or `WeightedTimeEdge`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> e1 = TimeEdge("A", "B", 1)
A->B at time 1

julia> edge_timestamp(e1)
1

julia> a = Edge(1, 2)
1->2

julia> edge_timestamp(a)
ERROR: type Edge has no field timestamp
[...]
```
"""
edge_timestamp(e::AbstractEdge) = e.timestamp

"""
    Edge(source, target)

Construct an Edge from a source node and a target node.
"""
struct Edge{V} <: AbstractEdge{V}
    source::V
    target::V
end

==(e1::Edge, e2::Edge) = (e1.source == e2.source && e1.target == e2.target)
eltype{V}(::Edge{V}) = V
eltype{V}(::Type{Edge{V}}) = V

"""
    TimeEdge(source, target, timestamp)

Construct a TimeEdge with source node `source`, target node `target` at time stamp `timestamp`.
"""
struct TimeEdge{V,T} <: AbstractEdge{V}
    source::V
    target::V
    timestamp::T
end

==(e1::TimeEdge, e2::TimeEdge) = (e1.source == e2.source &&
                                  e1.target == e2.target &&
                                  e1.timestamp == e2.timestamp)
eltype{V,T}(::TimeEdge{V,T}) = (V,T)
eltype{V,T}(::Type{TimeEdge{V,T}}) = (V,T)

"""
    WeightedTimeEdge(source, target, weight, timestamp)
    WeightedTimeEdge(source, target, timestamp)

Construct a WeightedTimeEdge. if `weight` is not given, set `weight = 1.0`.
"""
struct WeightedTimeEdge{V, T, W<:Real} <:AbstractEdge{V}
    source::V
    target::V
    weight::W
    timestamp::T
end
WeightedTimeEdge(source, target, timestamp) = WeightedTimeEdge(source, target, 1.0, timestamp)
==(e1::WeightedTimeEdge, e2::WeightedTimeEdge) =
                (e1.source == e2.source &&
                 e1.target == e2.target &&
                 e1.weight ≈ e2.weight &&
                 e1.timestamp == e2.timestamp)

eltype{V,T,W}(::WeightedTimeEdge{V,T,W}) = (V,T,W)
eltype{V,T,W}(::Type{WeightedTimeEdge{V,T,W}}) = (V,T,W)

"""
    edge_weight(e)

Returns the edge weight of a WeightedTimeEdge `e`.
"""
edge_weight(e::WeightedTimeEdge) = e.weight

"""
    edge_reverse(e)

Return the reverse of edge `e`, i.e., source is now target and target is now source.
"""
edge_reverse(e::Edge) = Edge(e.target, e.source)
edge_reverse(e::TimeEdge) = TimeEdge(e.target, e.source, e.timestamp)
edge_reverse(e::WeightedTimeEdge) =
      WeightedTimeEdge(e.target, e.source, e.weight, e.timestamp)

"""
    has_node(e, v)
    has_node(g, v)

In the first case, return `true` if `v` is a node of the edge `e`. In the second case, return `true` if graph `g` contains node `v` and `false` otherwise, where `v` can a node type object or a node key.
"""
has_node{V}(e::AbstractEdge{V}, v::V) = (v == source(e) || v == target(e))

"""
    has_active_node(g,v)
    has_active_node(g,key,t)

Return `true` if evolving graph `g` contains active node `v`
and `false` otherwise, where `v` is a TimeNode. Alternatively,
we can input node key `key` and timestamp `t` to check if `g`
contains the corresponding active node.
"""
function has_active_node end

"""
    AbstractPath

Abstract supertype for all path types.
"""
abstract type AbstractPath end

TimeOrTuple = Union{TimeNode,Tuple{Int,Int}}
"""
    TemporalPath()

Construct a TemporalPath, i.e., an array of TimeNode.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> p = TemporalPath()

julia> push!(p, TimeNode(1, "a", 2001))
 TimeNode(a, 2001)

julia> push!(p, TimeNode(1, "b", 2002))
TimeNode(a, 2001)->TimeNode(b, 2002)

julia> p
TimeNode(a, 2001)->TimeNode(b, 2002)

julia> append!(p,p)
TimeNode(a, 2001)->TimeNode(b, 2002)->TimeNode(a, 2001)->TimeNode(b, 2002)
```
"""
mutable struct TemporalPath{T<:TimeOrTuple} <: AbstractPath
    nodes::Vector{T}
end

TemporalPath() = TemporalPath(Union{TimeNode,Tuple{Int,Int}}[])
push!(t::TemporalPath, n::TimeOrTuple) = (push!(t.nodes, n);  t)
append!(t::TemporalPath, ns::Vector{TimeOrTuple}) = (append!(t.nodes, ns);  t)
append!(t::TemporalPath, t2::TemporalPath) = (append!(t, t2.nodes);  t)
deepcopy(t::TemporalPath) = TemporalPath(deepcopy(t.nodes))

length(p::TemporalPath) = length(p.nodes)
Base.start(P::TemporalPath) = 1
Base.next(p::TemporalPath, state) = (p.nodes[state], state+1)
Base.done(p::TemporalPath, state) = state > length(p.nodes)

==(p1::TemporalPath, p2::TemporalPath) = (p1.nodes == p2.nodes)

"""
    nodes(g)

Return the nodes of a graph `g`.
"""
function nodes end

"""
    num_nodes(g)
"""
function num_nodes end

"""
    active_nodes(g)

Return the active nodes of an evolving graph `g`.
"""
function active_nodes end

"""
    num_active_nodes(g)

Return the number of active nodes of evolving graph `g`.
"""
function num_active_nodes end

"""
    find_node(g, v)

Return node `v` if `v` is a node of graph `g`, otherwise return false. If `v` is a node key, return corresponding node.
"""
function find_node end

"""
    add_node!(g, v)

Add a node to graph `g`, where `v` can be either a node type object or a node key.
"""
function add_node! end

"""
    edges(g)
    edges(g, t)

Return all the edges of graph `g`. If timestamp `t` is given, returns all the edges at timestamp `t`.
"""
function edges end

"""
    num_edges(g)

Return the number of edges of graph `g`.
"""
function num_edges end

"""
    add_edge!(g, v1, v2, t; weight = 1.0)

Add an edge from `v1` to `v2` at timestamp `t` to evolving graph `g`. By default edge weight `weight` is equal to `1.0`.
"""
function add_edge! end

"""
    add_bunch_of_edges!(g, ebunch)

Add a bunch of edges to graph `g` where `ebunch` is an array of edges.
Each edge in `ebunch` is of form `(source, target, timestamp)` or
`(source, target, timestamp, weight)`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = EvolvingGraph()
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> add_bunch_of_edges!(g, [(1,2,2001), (2,3,2001), (3,1,2002), (2,3,2004)])
Directed EvolvingGraph 3 nodes, 4 static edges, 3 timestamps

julia> timestamps(g)
4-element Array{Int64,1}:
 2001
 2001
 2002
 2004
```
"""
function add_bunch_of_edges! end

"""
    is_directed(g)

Determine if a graph `g` is a directed graph.
"""
function is_directed end

"""
    timestamps(g)

Return the timestamps of graph `g`.
"""
function timestamps end

"""
    unique_timestamps(g)

Return the unique timestamps of graph `g`.
"""
function unique_timestamps end

"""
    num_timestamps(g)

Return the number of timestamps of graph `g`.
"""
function num_timestamps end

"""
    forward_neighbors(g, v)

Find the forward neighbors of a node `v` in graph `g`. If `g` is an evolving graph, we define the forward neighbors of a TimeNode `v` to be a collection of forward neighbors at time stamp `node_timestamp(v)` and the same node key at later time stamps.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = evolving_graph_from_arrays([1,1,2],[2,3,3], [1,2,3])
Directed EvolvingGraph 3 nodes, 3 static edges, 3 timestamps

julia> n = active_nodes(g)[1]
TimeNode(1, 1)

julia> forward_neighbors(g,n)
2-element Array{EvolvingGraphs.TimeNode{Int64,Int64},1}:
 TimeNode(2, 1)
 TimeNode(1, 2)
```

# References:

1. Jiahao Chen and Weijian Zhang, The Right Way to Search Evolving Graphs, Proceedings of IPDPS 2016.
"""
function forward_neighbors end

"""
    backward_neighbors(g, v)

Find the backward neighbors of a node `v` in graph `g`. If `g` is an evolving graph, we define the backward neighbors of a TimeNode `v` to be a collection of backward neighbors at time stamp `node_timestamp(v)` and the same node key at earlier time stamps.

```jldoctest
julia> using EvolvingGraphs

julia> g = evolving_graph_from_arrays([1,1,2],[2,3,3], [1,2,3])
Directed EvolvingGraph 3 nodes, 3 static edges, 3 timestamps

julia> n = active_nodes(g)[2]
TimeNode(2, 1)

julia> backward_neighbors(g,n)
1-element Array{EvolvingGraphs.TimeNode{Int64,Int64},1}:
 TimeNode(1, 1)
```

# References:

1. Jiahao Chen and Weijian Zhang, The Right Way to Search Evolving Graphs, Proceedings of IPDPS 2016.
"""
function backward_neighbors end


text:
# Copyright (c) 2015 Dahua Lin, Miles Lubin, Joey Huchette, Iain Dunning, and
#   contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using Gurobi
using Test

# When adding new tests, be mindful about creating new environments. Either
# re-use an existing environment in the module, or name the test function
# `test_MULTI_ENV_xxx` to trap the specific Gurobi error indicating that an
# environment could not be created.
const GRB_ENV = Gurobi.Env(output_flag = 0)

@testset "MathOptInterface Tests" begin
    for file in readdir("MOI")
        include(joinpath("MOI", file))
    end
end


text:
module Strings

using Compat, Mmap

# this would be handled by native Julia GC, but we'll roll our own for now
const PAGESIZE = @compat Int(@unix ? ccall(:jl_getpagesize, Clong, ()) : ccall(:jl_getallocationgranularity, Clong, ()))

type StringPool
    pool::Vector{Vector{UInt8}}
    ind::Int
    pos::Int
end

const POOL = StringPool(Any[Mmap.mmap(UInt8,PAGESIZE)],1,1)

function ensureroom!(n::Int)
    if POOL.pos + n < PAGESIZE
        # we have enough room to allocate `n` bytes
        return
    elseif n < PAGESIZE
        # we're hitting a page boundary
        push!(POOL.pool,Mmap.mmap(UInt8,PAGESIZE))
        POOL.ind += 1
        POOL.pos = 1
        return
    elseif n > PAGESIZE
        totalneededbytes = (div(n, PAGESIZE) + 1) * PAGESIZE
        push!(POOL.pool,Mmap.mmap(UInt8,totalneededbytes))
        POOL.ind += 1
        POOL.pos = 1
        return
    end
end

function storebytes!(s::Ptr{UInt8},n::Int,offset::Int=0)
    ensureroom!(n)
    # unsafe_copy!(dest::Ptr{T}, src::Ptr{T}, N)
    ptr = pointer(POOL.pool[POOL.ind])+UInt(POOL.pos)
    unsafe_copy!(ptr, s+UInt(offset), n)
    return ptr, n
end

storebytes!(s::Vector{UInt8},offset::Int=0) = storebytes!(pointer(s), length(s), offset)
storebytes!(s::Vector{UInt16}) = storebytes!(reinterpret(UInt8,s),length(s)*2)
storebytes!(s::Vector{UInt32}) = storebytes!(reinterpret(UInt8,s),length(s)*4)
storebytes!(s::Vector{Char})   = storebytes!(reinterpret(UInt8,s),length(s)*4)

# implementation files
include("encodings.jl")
include("string.jl")
include("ascii.jl")

end # module

macro s_str(s)
    Strings.String(s)
end


text:
using StatsBase
using Base.Test

##### Location

## geomean

@test_approx_eq geomean([1, 2, 3]) (6.0)^(1/3)
@test_approx_eq geomean(1:3) (6.0)^(1/3)
@test_approx_eq geomean([2, 8]) 4.0
@test_approx_eq geomean([4, 1, 1/32]) 0.5
@test geomean([1, 0, 2]) == 0.0

## harmmean

@test_approx_eq harmmean([1, 2, 3]) 3 / (1 + 1/2 + 1/3)
@test_approx_eq harmmean(1:3) 3 / (1 + 1/2 + 1/3)
@test_approx_eq harmmean([1, 2, 4]) 12 / 7

## genmean
@test_approx_eq genmean([1,1,2,3], 1) 7/4
@test_approx_eq genmean([1,4,2], -1) 12/7
@test_approx_eq genmean([1,1,2,3], 0) (6.0)^(1/4)
@test_approx_eq genmean([1.2,-0.5,0], 2) sqrt(169/300)
@test_approx_eq genmean([16/9,0.25,1.0], 1.5) (755/648)^(2/3)
# Test numerical stability for `p` close to 0 (genmean should be close to geometric mean).
@test_approx_eq_eps genmean([1,1,2,3], -1e-8) (6.0)^(1/4) 1e-8
# Test numerical stability for large `p` (genmean should be close to max).
@test_approx_eq_eps genmean([0.98,1.02], 1e4) 1.02 1e-4

## trimmean

@test_approx_eq trimmean([-100, 2, 3, 7, 200], 0.0) 22.4
@test_approx_eq trimmean([-100, 2, 3, 7, 200], 0.4) 4.0
@test_approx_eq trimmean([-100, 2, 3, 7, 200], 0.8) 3.0
@test_approx_eq trimmean([2, 3, -100, 200, 7], 0.4) 4.0
@test_approx_eq trimmean([2, 3, -100, 200, 7], 0.8) 3.0

## mode & modes

@test mode([1, 2, 3, 3, 2, 2, 1], 1:3) == 2
@test modes([1, 2, 3, 3, 2, 2, 1], 1:3) == [2]
@test modes([1, 3, 2, 3, 3, 2, 2, 1], 1:3) == [2, 3]

@test mode([1, 2, 3, 3, 2, 2, 1]) == 2
@test modes([1, 2, 3, 3, 2, 2, 1]) == [2]
@test sort(modes([1, 3, 2, 3, 3, 2, 2, 1])) == [2, 3]

## zscores

@test zscore([-3:3;], 1.5, 0.5) == [-9.0:2.0:3.0;]

a = [3 4 5 6; 7 8 1 2; 6 9 3 0]
z1 = [4. 6. 8. 10.; 5. 6. -1. 0.; 1.5 3.0 0.0 -1.5]
z2 = [8. 2. 3. 1.; 24. 10. -1. -1.; 20. 12. 1. -2.]

@test_approx_eq zscore(a, [1, 2, 3], [0.5, 1.0, 2.0]) z1
@test_approx_eq zscore(a, [1 3 2 4], [0.25 0.5 1.0 2.0]) z2

@test zscore!(collect(-3.0:3.0), 1.5, 0.5) == [-9.0:2.0:3.0;]
@test_approx_eq zscore!(float(a), [1, 2, 3], [0.5, 1.0, 2.0]) z1
@test_approx_eq zscore!(float(a), [1 3 2 4], [0.25 0.5 1.0 2.0]) z2

@test zscore!(zeros(7), [-3:3;], 1.5, 0.5) == [-9.0:2.0:3.0;]
@test_approx_eq zscore!(zeros(size(a)), a, [1, 2, 3], [0.5, 1.0, 2.0]) z1
@test_approx_eq zscore!(zeros(size(a)), a, [1 3 2 4], [0.25 0.5 1.0 2.0]) z2

@test_approx_eq zscore(a) zscore(a, mean(a), std(a))
@test_approx_eq zscore(a, 1) zscore(a, mean(a,1), std(a,1))
@test_approx_eq zscore(a, 2) zscore(a, mean(a,2), std(a,2))

###### quantile & friends

@test_approx_eq quantile(1:5) [1:5;]
@test_approx_eq nquantile(1:5, 2) [1, 3, 5]
@test_approx_eq nquantile(1:5, 4) [1:5;]

@test_approx_eq percentile([1:5;], 25) 2.0
@test_approx_eq percentile([1:5;], [25, 50, 75]) [2.0, 3.0, 4.0]

##### Dispersion

@test span([3, 4, 5, 6, 2]) == (2:6)

@test_approx_eq variation([1:5;]) 0.527046276694730

@test_approx_eq sem([1:5;]) 0.707106781186548

@test StatsBase.mad!([1:5;], 3) ≈ 1.4826022185056018
@test mad(1:5) ≈ 1.4826022185056018
@test StatsBase.mad!([1:5;], constant=1.0) ≈ 1.0
@test StatsBase.mad!([1:5;], 3, constant=1.0) ≈ 1.0
@test_throws ArgumentError mad(Int[])

# Issue 197
@test mad(1:2) ≈ 0.7413011092528009

@test_approx_eq iqr(1:5) 2.0

##### entropy

@test_approx_eq entropy([0.5, 0.5]) 0.6931471805599453
@test_approx_eq entropy([0.2, 0.3, 0.5]) 1.0296530140645737

@test_approx_eq entropy([0.5, 0.5],2) 1.0
@test_approx_eq entropy([0.2, 0.3, 0.5], 2) 1.4854752972273344

##### Renyi entropies
# Generate a random probability distribution
nindiv = 50
dist = rand(nindiv)
dist /= sum(dist)

# Check Shannon entropy against Renyi entropy of order 1
@test_approx_eq entropy(dist) renyientropy(dist, 1)
@test_approx_eq renyientropy(dist, 1) renyientropy(dist, 1.0)

# Check Renyi entropy of order 0 is the natural log of the count of non-zeros
@test_approx_eq renyientropy(dist, 0) log(mapreduce(x -> x > 0 ? 1 : 0, +, dist))

# And is therefore not affected by the addition of non-zeros
zdist = dist
zdist = append!(dist, zeros(rand(50)))
@test_approx_eq renyientropy(dist, 0) renyientropy(zdist, 0)

# Indeed no Renyi entropy should be
loworder = rand() # Between 0 and 1
@test_approx_eq renyientropy(dist, loworder) renyientropy(zdist, loworder)
order = rand() * 49 + 1 # And over 1
@test_approx_eq renyientropy(dist, order) renyientropy(zdist, order)

# Renyi entropy of order infinity is -log(maximum(dist))
@test_approx_eq renyientropy(dist, Inf) -log(maximum(dist))

# And all Renyi entropies of uniform probabilities should be log n, where n is the length
udist = ones(nindiv) / nindiv
@test_approx_eq renyientropy(udist, order) log(nindiv)

# And test generalised probability distributions (sum(p) != 1)
scale = rand()
@test_approx_eq renyientropy(udist * scale, 0) renyientropy(udist, 0) - log(scale)
@test_approx_eq renyientropy(udist * scale, 1) renyientropy(udist, 1) - log(scale)
@test_approx_eq renyientropy(udist * scale, Inf) renyientropy(udist, Inf) - log(scale)
@test_approx_eq renyientropy(udist * scale, order) renyientropy(udist, order) - log(scale)

##### Cross entropy
@test_approx_eq crossentropy([0.2, 0.3, 0.5], [0.3, 0.4, 0.3]) 1.1176681825904018
@test_approx_eq crossentropy([0.2, 0.3, 0.5], [0.3, 0.4, 0.3], 2) 1.6124543443825532

##### KL divergence
@test_approx_eq kldivergence([0.2, 0.3, 0.5], [0.3, 0.4, 0.3]) 0.08801516852582819
@test_approx_eq kldivergence([0.2, 0.3, 0.5], [0.3, 0.4, 0.3], 2) 0.12697904715521868

##### summarystats

s = summarystats(1:5)
@test isa(s, StatsBase.SummaryStats)
@test s.min == 1.0
@test s.max == 5.0
@test_approx_eq s.mean 3.0
@test_approx_eq s.median 3.0
@test_approx_eq s.q25 2.0
@test_approx_eq s.q75 4.0


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Chemfiles"
version = v"0.10.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/chemfiles/chemfiles/archive/$version.tar.gz",
                  "a1bbc26d8e8c33f52a157df4aece0c07e126a986f29231c4fdcb9a55292dc5b0"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/chemfiles-*/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libchemfiles", :libchemfiles)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
@testset "utility.jl" begin
    p1 = rad2deg(Point(0.5, 1.0))
    p2 = rad2deg(Point(1.0, 0.5))
    p3 = rad2deg(Point(0.6, 1.0))  # above p1

    rs = RouteSection(p1,p2)

    @test (p1-p2).ϕ ≈ -28.6478897565 atol = 0.001

    @test Navigation.normalize(-270.0, -180.0, 180.0) ≈ 90.0 atol = 0.1
    @test Navigation.normalize(181.0, -180.0, 180.0) ≈ -179.0 atol = 0.1
    @test deg2rad(p1).ϕ ≈ 0.5 atol = 0.0001
    @test deg2rad(p2).λ ≈ 0.5 atol = 0.0001

    p3 = Point(80.0, 190.0)
    @test Navigation.normalize(p3).ϕ ≈ rad2deg(1.3962634016) atol = 0.0001
    @test Navigation.normalize(p3).λ ≈ rad2deg(-2.96705972839) atol = 0.0001
    @test Navigation.normalize(p3).λ ≈ -170.000 atol = 0.01

end