text:
using Dates

using Pumas, CairoMakie, PumasPlots, GlobalSensitivity

model = @model begin
    @param begin
        Fup ∈ RealDomain(init = 0.42)
        fumic ∈ RealDomain(init = 0.711)
        WEIGHT ∈ RealDomain(init = 73)
        MPPGL ∈ RealDomain(init = 30.3)
        MPPGI ∈ RealDomain(init = 0)
        C_OUTPUT ∈ RealDomain(init = 6.5)
        VmaxH ∈ RealDomain(init = 40)
        VmaxG ∈ RealDomain(init = 40)
        KmH ∈ RealDomain(init = 9.3)
        KmG ∈ RealDomain(init = 9.3)
        bp ∈ RealDomain(init = 1)
        kpad ∈ RealDomain(init = 9.89)
        kpbo ∈ RealDomain(init = 7.91)
        kpbr ∈ RealDomain(init = 7.35)
        kpgu ∈ RealDomain(init = 5.82)
        kphe ∈ RealDomain(init = 1.95)
        kpki ∈ RealDomain(init = 2.9)
        kpli ∈ RealDomain(init = 4.66)
        kplu ∈ RealDomain(init = 0.83)
        kpmu ∈ RealDomain(init = 2.94)
        kpsp ∈ RealDomain(init = 2.96)
        kpre ∈ RealDomain(init = 4)
        MW ∈ RealDomain(init = 349.317)
        logP ∈ RealDomain(init = 2.56)
        s_lumen ∈ RealDomain(init = 0.39*1000)
        L ∈ RealDomain(init = 280)
        d ∈ RealDomain(init = 2.5)
        PF ∈ RealDomain(init = 1.57)
        VF ∈ RealDomain(init = 6.5)
        MF ∈ RealDomain(init = 13)
        ITT ∈ RealDomain(init = 3.32)
        A ∈ RealDomain(init = 7440)
        B ∈ RealDomain(init = 1e7)
        alpha ∈ RealDomain(init = 0.6)
        beta ∈ RealDomain(init = 4.395)
        fabs ∈ RealDomain(init = 1)
        fdis ∈ RealDomain(init = 1)
        fperm ∈ RealDomain(init = 1)
        vad ∈ RealDomain(init = 18.2)
        vbo ∈ RealDomain(init =10.5)
        vbr ∈ RealDomain(init =1.45)
        vguWall ∈ RealDomain(init =0.65)
        vgulumen ∈ RealDomain(init =0.35)
        vhe ∈ RealDomain(init =0.33)
        vki ∈ RealDomain(init =0.31)
        vli ∈ RealDomain(init =1.8)
        vlu ∈ RealDomain(init =0.5)
        vmu ∈ RealDomain(init =29)
        vsp ∈ RealDomain(init =0.15)
        vbl ∈ RealDomain(init =5.6)
        FQad ∈ RealDomain(lower = 0.0, init = 0.05, upper = 1.0) #add bounds to parameters for estimation
        FQbo ∈ RealDomain(lower = 0.0, init = 0.05, upper = 1.0)
        FQbr ∈ RealDomain(lower = 0.0, init = 0.12, upper = 1.0)
        FQgu ∈ RealDomain(lower = 0.0, init = 0.16, upper = 1.0)
        FQhe ∈ RealDomain(lower = 0.0, init = 0.04, upper = 1.0)
        FQki ∈ RealDomain(lower = 0.0, init = 0.19, upper = 1.0)
        FQli ∈ RealDomain(lower = 0.0, init = 0.255, upper = 1.0)
        FQmu ∈ RealDomain(lower = 0.0, init = 0.17, upper = 1.0)
        FQsp ∈ RealDomain(lower = 0.0, init = 0.03, upper = 1.0)
    end
    @pre begin
        Vgu = vguWall + vgulumen
        Vve = 0.705*vbl
        Var = 0.295*vbl
        Vre = WEIGHT - (vli+vki+vsp+vhe+vlu+vbo+vbr+vmu+vad+vguWall+vbl)
        CO = C_OUTPUT*60
        Qad = FQad*CO
        Qbo = FQbo*CO
        Qbr = FQbr*CO
        Qgu = FQgu*CO
        Qhe = FQhe*CO
        Qki = FQki*CO
        Qli = FQli*CO
        Qmu = FQmu*CO
        Qsp = FQsp*CO
        Qha = Qli - (Qgu+Qsp)
        Qtot = Qli+Qki+Qbo+Qhe+Qmu+Qad+Qbr
        Qre = CO - Qtot
        Qlu = CO
        Vgulumen = vgulumen
        S_lumen = s_lumen
        VguWall = vguWall
        Kpgu = kpgu
        BP = bp
        Vad = vad
        Kpad = kpad
        Vbr = vbr
        Kpbr = kpbr
        Vhe = vhe
        Kphe = kphe
        Vki = vki
        Kpki = kpki
        fup = Fup
        Vsp = vsp
        Kpsp = kpsp
        Vli = vli
        Kpli = kpli
        Vlu = vlu
        Kplu = kplu
        Kpmu = kpmu
        Kpre = kpre
        Vmu = vmu
        Vbl = vbl
        Vbo = vbo
        Kpbo = kpbo
        SA_abs = pi*L*d*PF*VF*MF*1e-4
        SA_basal = pi*L*d*PF*VF*1e-4
        MA = 10^logP
        MW_eff = MW - (3*17)
        Peff = fperm*A*(((MW_eff^(-alpha-beta))*MA)/((MW_eff^(-alpha)) + B*(MW_eff^(-beta))*MA) * 1e-2 * 3600)
        kd = fdis*Peff*SA_abs*1000/vgulumen
        ka = fabs*Peff*SA_basal*1000/VguWall
        kt = 1/ITT
        scale_factor_H = MPPGL*Vli*1000
        scale_factor_G = MPPGI*VguWall*1000
        CLintHep = ((VmaxH/KmH)*scale_factor_H*60*1e-6)/fumic
        CLintGut = ((VmaxG/KmG)*scale_factor_G*60*1e-6)/fumic
        #CLintHep = CLintHep/fumic
        #CLintGut = CLintGut/fumic
        CLrenal = 0.096
        f = 1
    end
    @dynamics begin
        GUTLUMEN' = -kd*Vgulumen*(f*(GUTLUMEN/Vgulumen) + (1-f)*S_lumen) -
            kt*GUTLUMEN
        GUTWALL' = kd*Vgulumen*(f*(GUTLUMEN/Vgulumen) + (1-f)*S_lumen) -
            ka*GUTWALL - CLintGut*(GUTWALL/VguWall)
        GUT' = ka*GUTWALL + Qgu*((ART/Var) - (GUT/VguWall)/(Kpgu/BP))
        ADIPOSE' = Qad*((ART/Var) - (ADIPOSE/Vad)/(Kpad/BP))
        BRAIN' = Qbr*((ART/Var) - (BRAIN/Vbr)/(Kpbr/BP))
        HEART' = Qhe*((ART/Var) - (HEART/Vhe)/(Kphe/BP))
        KIDNEY' = Qki*((ART/Var) - (KIDNEY/Vki)/(Kpki/BP)) -
            CLrenal*(((KIDNEY/Vki)*fup)/(Kpki/BP))
        LIVER' = Qgu*((GUT/VguWall)/(Kpgu/BP)) + Qsp*((SPLEEN/Vsp)/(Kpsp/BP)) +
            Qha*(ART/Var) - Qli*((LIVER/Vli)/(Kpli/BP)) -
            CLintHep*(((LIVER/Vli)*fup)/(Kpli/BP))
        LUNG' = Qlu*((VEN/Vve) - (LUNG/Vlu)/(Kplu/BP))
        MUSCLE' = Qmu*((ART/Var) - (MUSCLE/Vmu)/(Kpmu/BP))
        SPLEEN' = Qsp*((ART/Var) - (SPLEEN/Vsp)/(Kpsp/BP))
        BONE' = Qbo*((ART/Var) - (BONE/Vbo)/(Kpbo/BP))
        REST' = Qre*((ART/Var) - (REST/Vre)/(Kpre/BP))
        VEN' = Qad*((ADIPOSE/Vad)/(Kpad/BP)) + Qbr*((BRAIN/Vbr)/(Kpbr/BP)) +
            Qhe*((HEART/Vhe)/(Kphe/BP)) + Qki*((KIDNEY/Vki)/(Kpki/BP)) +
            Qli*((LIVER/Vli)/(Kpli/BP)) + Qmu*((MUSCLE/Vmu)/(Kpmu/BP)) +
            Qbo*((BONE/Vbo)/(Kpbo/BP)) + Qre*((REST/Vre)/(Kpre/BP)) -
            Qlu*(VEN/Vve)
        ART' = Qlu*((LUNG/Vlu)/(Kplu/BP) - (ART/Var))
    end
    @derived begin
        Cvenn = VEN./Vve
        cp ~ @. Normal(Cvenn, 0.1) #for estimation
    end
end

regimen_s = DosageRegimen(200, time=0, addl=13, ii=12, cmt=1, ss=1)
sub_s = Subject(id=1, events=regimen_s)

p = (Fup = 0.42, fumic = 0.711, WEIGHT = 73, MPPGL = 30.3, MPPGI = 0,
    C_OUTPUT = 6.5, VmaxH = 40, VmaxG = 40, KmH = 9.3, KmG = 9.3, bp = 1,
    kpad = 9.89, kpbo = 7.91, kpbr = 7.35, kpgu = 5.82, kphe = 1.95, kpki = 2.9,
    kpli = 4.66, kplu = 0.83, kpmu = 2.94, kpsp = 2.96, kpre = 4, MW = 349.317,
    logP = 2.56, s_lumen = 0.39*1000, L = 280, d = 2.5, PF = 1.57, VF = 6.5,
    MF = 13, ITT = 3.32, A = 7440, B = 1e7, alpha = 0.6, beta = 4.395, fabs = 1,
    fdis = 1, fperm = 1, vad = 18.2, vbo = 10.5, vbr = 1.45, vguWall = 0.65,
    vgulumen = 0.35, vhe = 0.33, vki = 0.31, vli = 1.8, vlu = 0.5, vmu = 29,
    vsp = 0.15, vbl = 5.6, FQad = 0.05, FQbo = 0.05, FQbr = 0.12, FQgu = 0.16,
    FQhe = 0.04, FQki = 0.19, FQli = 0.255, FQmu = 0.17, FQsp = 0.03)

simdata = simobs(model, [sub_s], p)
sim_plot(model, simdata, observations=[:Cvenn])

data = read_pumas(DataFrame(simdata), observations = [:cp])
ft = fit(model, data, p, Pumas.NaivePooled(),
    constantcoef = (
        Fup = 0.42, fumic = 0.711, WEIGHT = 73, MPPGL = 30.3, MPPGI = 0,
        C_OUTPUT = 6.5, VmaxH = 40, VmaxG = 40, KmH = 9.3, KmG = 9.3, bp = 1,
        kpad = 9.89, kpbo = 7.91, kpbr = 7.35, kpgu = 5.82, kphe = 1.95,
        kpki = 2.9, kpli = 4.66, kplu = 0.83, kpmu = 2.94, kpsp = 2.96,
        kpre = 4, MW = 349.317, logP = 2.56, s_lumen = 0.39*1000, L = 280,
        d = 2.5, PF = 1.57, VF = 6.5, MF = 13, ITT = 3.32, A = 7440, B = 1e7,
        alpha = 0.6, beta = 4.395, fabs = 1, fdis = 1, fperm = 1, vad = 18.2,
        vbo = 10.5, vbr = 1.45, vguWall = 0.65, vgulumen = 0.35, vhe = 0.33,
        vki = 0.31, vli = 1.8, vlu = 0.5, vmu = 29, vsp = 0.15, vbl = 5.6),
    ensemblealg=EnsembleThreads())

model = @model begin
    @param begin
        Fup ∈ RealDomain(init = 0.42)
        fumic ∈ RealDomain(init = 0.711)
        WEIGHT ∈ RealDomain(init = 73)
        MPPGL ∈ RealDomain(init = 30.3)
        MPPGI ∈ RealDomain(init = 0)
        C_OUTPUT ∈ RealDomain(init = 6.5)
        VmaxH ∈ RealDomain(init = 40)
        VmaxG ∈ RealDomain(init = 40)
        KmH ∈ RealDomain(init = 9.3)
        KmG ∈ RealDomain(init = 9.3)
        bp ∈ RealDomain(init = 1)
        kpad ∈ RealDomain(init = 9.89)
        kpbo ∈ RealDomain(init = 7.91)
        kpbr ∈ RealDomain(init = 7.35)
        kpgu ∈ RealDomain(init = 5.82)
        kphe ∈ RealDomain(init = 1.95)
        kpki ∈ RealDomain(init = 2.9)
        kpli ∈ RealDomain(init = 4.66)
        kplu ∈ RealDomain(init = 0.83)
        kpmu ∈ RealDomain(init = 2.94)
        kpsp ∈ RealDomain(init = 2.96)
        kpre ∈ RealDomain(init = 4)
        MW ∈ RealDomain(init = 349.317)
        logP ∈ RealDomain(init = 2.56)
        s_lumen ∈ RealDomain(init = 0.39*1000)
        L ∈ RealDomain(init = 280)
        d ∈ RealDomain(init = 2.5)
        PF ∈ RealDomain(init = 1.57)
        VF ∈ RealDomain(init = 6.5)
        MF ∈ RealDomain(init = 13)
        ITT ∈ RealDomain(init = 3.32)
        A ∈ RealDomain(init = 7440)
        B ∈ RealDomain(init = 1e7)
        alpha ∈ RealDomain(init = 0.6)
        beta ∈ RealDomain(init = 4.395)
        fabs ∈ RealDomain(init = 1)
        fdis ∈ RealDomain(init = 1)
        fperm ∈ RealDomain(init = 1)
        vad ∈ RealDomain(init = 18.2)
        vbo ∈ RealDomain(init =10.5)
        vbr ∈ RealDomain(init =1.45)
        vguWall ∈ RealDomain(init =0.65)
        vgulumen ∈ RealDomain(init =0.35)
        vhe ∈ RealDomain(init =0.33)
        vki ∈ RealDomain(init =0.31)
        vli ∈ RealDomain(init =1.8)
        vlu ∈ RealDomain(init =0.5)
        vmu ∈ RealDomain(init =29)
        vsp ∈ RealDomain(init =0.15)
        vbl ∈ RealDomain(init =5.6)
        FQad ∈ RealDomain(lower = 0.0, init = 0.05, upper = 1.0)
        FQbo ∈ RealDomain(lower = 0.0, init = 0.05, upper = 1.0)
        FQbr ∈ RealDomain(lower = 0.0, init = 0.12, upper = 1.0)
        FQgu ∈ RealDomain(lower = 0.0, init = 0.16, upper = 1.0)
        FQhe ∈ RealDomain(lower = 0.0, init = 0.04, upper = 1.0)
        FQki ∈ RealDomain(lower = 0.0, init = 0.19, upper = 1.0)
        FQli ∈ RealDomain(lower = 0.0, init = 0.255, upper = 1.0)
        FQmu ∈ RealDomain(lower = 0.0, init = 0.17, upper = 1.0)
        FQsp ∈ RealDomain(lower = 0.0, init = 0.03, upper = 1.0)
    end
    @pre begin
        Vgu = vguWall + vgulumen
        Vve = 0.705*vbl
        Var = 0.295*vbl
        Vre = WEIGHT - (vli+vki+vsp+vhe+vlu+vbo+vbr+vmu+vad+vguWall+vbl)
        CO = C_OUTPUT*60
        Qad = FQad*CO
        Qbo = FQbo*CO
        Qbr = FQbr*CO
        Qgu = FQgu*CO
        Qhe = FQhe*CO
        Qki = FQki*CO
        Qli = FQli*CO
        Qmu = FQmu*CO
        Qsp = FQsp*CO
        Qha = Qli - (Qgu+Qsp)
        Qtot = Qli+Qki+Qbo+Qhe+Qmu+Qad+Qbr
        Qre = CO - Qtot
        Qlu = CO
        Vgulumen = vgulumen
        S_lumen = s_lumen
        VguWall = vguWall
        Kpgu = kpgu
        BP = bp
        Vad = vad
        Kpad = kpad
        Vbr = vbr
        Kpbr = kpbr
        Vhe = vhe
        Kphe = kphe
        Vki = vki
        Kpki = kpki
        fup = Fup
        Vsp = vsp
        Kpsp = kpsp
        Vli = vli
        Kpli = kpli
        Vlu = vlu
        Kplu = kplu
        Kpmu = kpmu
        Kpre = kpre
        Vmu = vmu
        Vbl = vbl
        Vbo = vbo
        Kpbo = kpbo
        SA_abs = pi*L*d*PF*VF*MF*1e-4
        SA_basal = pi*L*d*PF*VF*1e-4
        MA = 10^logP
        MW_eff = MW - (3*17)
        Peff = fperm*A*(((MW_eff^(-alpha-beta))*MA)/((MW_eff^(-alpha)) + B*(MW_eff^(-beta))*MA) * 1e-2 * 3600)
        kd = fdis*Peff*SA_abs*1000/vgulumen
        ka = fabs*Peff*SA_basal*1000/VguWall
        kt = 1/ITT
        scale_factor_H = MPPGL*Vli*1000
        scale_factor_G = MPPGI*VguWall*1000
        CLintHep = ((VmaxH/KmH)*scale_factor_H*60*1e-6)/fumic
        CLintGut = ((VmaxG/KmG)*scale_factor_G*60*1e-6)/fumic
        #CLintHep = CLintHep/fumic
        #CLintGut = CLintGut/fumic
        CLrenal = 0.096
        f = 1
    end
    @dynamics begin
        GUTLUMEN' = -kd*Vgulumen*(f*(GUTLUMEN/Vgulumen) + (1-f)*S_lumen) -
            kt*GUTLUMEN
        GUTWALL' = kd*Vgulumen*(f*(GUTLUMEN/Vgulumen) + (1-f)*S_lumen) -
            ka*GUTWALL - CLintGut*(GUTWALL/VguWall)
        GUT' = ka*GUTWALL + Qgu*((ART/Var) - (GUT/VguWall)/(Kpgu/BP))
        ADIPOSE' = Qad*((ART/Var) - (ADIPOSE/Vad)/(Kpad/BP))
        BRAIN' = Qbr*((ART/Var) - (BRAIN/Vbr)/(Kpbr/BP))
        HEART' = Qhe*((ART/Var) - (HEART/Vhe)/(Kphe/BP))
        KIDNEY' = Qki*((ART/Var) - (KIDNEY/Vki)/(Kpki/BP)) -
            CLrenal*(((KIDNEY/Vki)*fup)/(Kpki/BP))
        LIVER' = Qgu*((GUT/VguWall)/(Kpgu/BP)) + Qsp*((SPLEEN/Vsp)/(Kpsp/BP)) +
            Qha*(ART/Var) - Qli*((LIVER/Vli)/(Kpli/BP)) -
            CLintHep*(((LIVER/Vli)*fup)/(Kpli/BP))
        LUNG' = Qlu*((VEN/Vve) - (LUNG/Vlu)/(Kplu/BP))
        MUSCLE' = Qmu*((ART/Var) - (MUSCLE/Vmu)/(Kpmu/BP))
        SPLEEN' = Qsp*((ART/Var) - (SPLEEN/Vsp)/(Kpsp/BP))
        BONE' = Qbo*((ART/Var) - (BONE/Vbo)/(Kpbo/BP))
        REST' = Qre*((ART/Var) - (REST/Vre)/(Kpre/BP))
        VEN' = Qad*((ADIPOSE/Vad)/(Kpad/BP)) + Qbr*((BRAIN/Vbr)/(Kpbr/BP)) +
            Qhe*((HEART/Vhe)/(Kphe/BP)) + Qki*((KIDNEY/Vki)/(Kpki/BP)) +
            Qli*((LIVER/Vli)/(Kpli/BP)) + Qmu*((MUSCLE/Vmu)/(Kpmu/BP)) +
            Qbo*((BONE/Vbo)/(Kpbo/BP)) + Qre*((REST/Vre)/(Kpre/BP)) -
            Qlu*(VEN/Vve)
        ART' = Qlu*((LUNG/Vlu)/(Kplu/BP) - (ART/Var))
    end
    @derived begin
        Cvenn = VEN./Vve
        #capturing NCA metrics for evaluations
        nca := @nca Cvenn
        auc =  last(NCA.auc(nca))
        cmax = last(NCA.cmax(nca))
    end
end

p_range_low = (fperm=1/3, s_lumen=390/3, ITT = 3.32/3, MPPGI=1.44/3, )

p_range_high = (fperm=1*3, s_lumen=390*3, ITT = 3.32*3, MPPGI=1.44*3, )

regimen_s = DosageRegimen(200, time=0, addl=13, ii=12, cmt=1, ss=1, route = Pumas.NCA.IVInfusion)
sub_s = Subject(id=1, events=regimen_s)
sobol_ = Pumas.gsa(model, sub_s, p, GlobalSensitivity.Sobol(), [:cmax,:auc], p_range_low,p_range_high, N=1000, obstimes=0.0:1.0:30.0)

keys_ = keys(p_range_low)
cmax_s1 = [sobol_.first_order[1,:][key] for key in keys_]
cmax_st = [sobol_.total_order[1,:][key] for key in keys_]

fig = Figure(resolution = (1200, 800))
plot_cmax_s1 = scatter(fig[1,1], 1:4, cmax_s1, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "First Order", title="Cmax"))
plot_cmax_st = scatter(fig[1,2], 1:4, cmax_st, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "Total Order"), marker=:utriangle)

auc_s1 = [sobol_.first_order[2,:][key] for key in keys_]
auc_st = [sobol_.total_order[2,:][key] for key in keys_]

plot_auc_s1 = scatter(fig[2,1], 1:4, auc_s1, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "First Order", title="AUC"))
plot_auc_st = scatter(fig[2,2], 1:4, auc_st, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "Total Order"), marker=:utriangle)
display(fig)

eFAST_ = Pumas.gsa(model, sub_s, p, GlobalSensitivity.eFAST(), [:cmax,:auc], p_range_low, p_range_high, n=1000, obstimes=0.0:1.0:30.0)

keys_ = keys(p_range_low)
cmax_s1 = [eFAST_.first_order[1,:][key] for key in keys_]
cmax_st = [eFAST_.total_order[1,:][key] for key in keys_]

fig = Figure(resolution = (1200,800))
plot_cmax_s1 = scatter(fig[1,1], 1:4, cmax_s1, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "First Order", title="Cmax"))
plot_cmax_st = scatter(fig[1,2], 1:4, cmax_st, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "Total Order"), marker=:utriangle)

auc_s1 = [eFAST_.first_order[2,:][key] for key in keys_]
auc_st = [eFAST_.total_order[2,:][key] for key in keys_]

plot_auc_s1 = scatter(fig[2,1], 1:4, auc_s1, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "First Order", title="AUC"))
plot_auc_st = scatter(fig[2,2], 1:4, auc_st, axis = (yticks = 0:1, xticks = (1:4, [string.(keys_)...]), label = "Total Order"), marker=:utriangle)
display(fig)


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.Execute
using QuantumCircuits.Execute: _0,  _1, qjacobian
using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Circuit: toQiskit
using QuantumCircuits.QCircuits.Circuit

using Zygote
using Zygote: jacobian

@test ket"0" == _0
@test ket"1" == _1

@test ket"00" == kron(_0, _0)
@test ket"01" == kron(_0, _1)
@test ket"10" == kron(_1, _0)
@test ket"11" == kron(_1, _1)

@test ket"0101" == kron(_0, kron(_1, kron(_0, _1)))

const backend = QuantumSimulator()
const qiskitBackend = QiskitQuantum()

#############################################################################

err(x, y) = sum(abs.(x - y))

function doTests(qc, expected, measuere_mat, checkQiskit=true)
    @test err(execute(backend, qc), expected) < 1e-5
    @test qc.measures_matrix == measuere_mat
    if checkQiskit
        @test err(execute(qiskitBackend, toQiskit(qc)), expected) < 0.05
    end
end

circuits = QCircuit[]

qc = QCircuit(2)
qc.h(0)
qc.x(1)
qc.measure([0, 1], [0, 1])
doTests(qc, [0.0, 0.0, 0.5, 0.5], eye(4))
push!(circuits, qc)

###

qc = QCircuit(2)
qc.h(0)
qc.x(1)
qc.measure(0, 0)
doTests(qc, [0.5, 0.5], [1 0 1 0; 0 1 0 1], false)

###

qr = QuantumRegister(2)
cr = ClassicalRegister(1)
qc = QCircuit(qr, cr)
qc.h(0)
qc.x(1)
qc.measure(0, 0)
doTests(qc, [0.5, 0.5], [1 0 1 0; 0 1 0 1])
push!(circuits, qc)

###

qc = QCircuit(2)
qc.h(0)
qc.x(1)
qc.measure(1, 0)
doTests(qc, [0.0, 1.0], [1 1 0 0; 0 0 1 1], false)

###

qc = QCircuit(2)
qc.h(0)
qc.x(1)
qc.measure(1, 1)
doTests(qc, [0.0, 1.0], [1 1 0 0; 0 0 1 1], false)

###

qr = QuantumRegister(2)
cr = ClassicalRegister(1)
qc = QCircuit(qr, cr)
qc.h(0)
qc.x(1)
qc.measure(1, 0)
doTests(qc, [0.0, 1.0], [1 1 0 0; 0 0 1 1])
push!(circuits, qc)

###

qr = QuantumRegister(3)
cr = ClassicalRegister(2)
qc = QCircuit(qr, cr)
qc.h(0)
qc.x(1)
qc.measure([0, 2], [0, 1])
expected = [0.5, 0.5, 0.0, 0.0]
@test err(execute(backend, qc), expected) < 1e-5
@test err(execute(qiskitBackend, toQiskit(qc)), expected) < 0.05
push!(circuits, qc)

###

qr = QuantumRegister(3)
cr = ClassicalRegister(2)
qc = QCircuit(qr, cr)
qc.h(0)
qc.x(1)
qc.measure([0, 1], [0, 1])
expected = [0.0, 0.0, 0.5, 0.5]
@test err(execute(backend, qc), expected) < 1e-5
@test err(execute(qiskitBackend, toQiskit(qc)), expected) < 0.05
push!(circuits, qc)

#############################################################################
res_sim = execute(backend, circuits)
res_dev = execute(qiskitBackend, [toQiskit(qc) for qc in circuits], debug=false)

@test sum([err(s, d) for (s, d) in zip(res_sim, res_dev)]) < 0.5

#############################################################################

qr = QuantumRegister(3)
cr = ClassicalRegister(2)
qc = QCircuit(qr, cr)
qc.rx(0)
qc.rx(1)
qc.rx(2)
qc.measure([0, 1], [0, 1])

pa = [1.0, 1.0, 1.0]
execute(backend, qc, pa)

loss(params::Vector{Float64}) = loss_expected_zero_state(execute(backend, qc, params))
loss(pa)
loss'(pa)

#############################################################################

@test loss_expected_zero_state(ket"0") < 1e-8
@test loss_expected_zero_state(ket"00") < 1e-8
@test loss_expected_zero_state(ket"000") < 1e-8
@test loss_expected_zero_state(ket"0000") < 1e-8

val1 = loss_expected_zero_state([1/2 1/6 1/6 1/6])
val2 = loss_expected_zero_state([1/2 1/2 0.0 0.0])
@test val1 + 0.1 < val2

val1 = loss_expected_zero_state([0.0 1/3 1/3 1/3])
val2 = loss_expected_zero_state([0.0 1.0 0.0 0.0])
@test val1 + 0.1 < val2

#############################################################################

qr = QuantumRegister(3)
cr = ClassicalRegister(3)
qc = QCircuit(qr, cr)
qc.rx(0)
qc.rx(1)
qc.rx(2)
qc.measure([0, 1, 2], [0, 1, 2])

pa = [1.0, 0.2, 0.6]

# Simulator
executeSim(params::Vector{Float64}) = execute(backend, qc, params)
lossAD(params::Vector{Float64}) = loss_expected_zero_state(executeSim(params))

res_sim = executeSim(pa)
res_sim_los = lossAD(pa)
res_sim_der = real.(lossAD'(pa))
res_sim_jak = jacobian(executeSim, pa)[1]

res_dev, res_dev_jak = qjacobian(qiskitBackend, qc, pa)
@test unitary_error(res_dev, res_sim) < 0.001
@test unitary_error(res_dev_jak, res_sim_jak) < 0.001

fn(params::Vector) = qderivative(qiskitBackend, qc, loss_expected_zero_state, params)
res_dev_los, res_dev_der = fn(pa)

@test unitary_error(res_dev_los, res_sim_los) < 0.001
@test unitary_error(res_dev_der, res_sim_der) < 0.001

#############################################################################

qr = QuantumRegister(2)
qc = QCircuit(qr)
qc.h(0)
qc.x(1)
@test err(execute(backend, qc), [0.0, 0.0, 0.5, 0.5]) < 1e-5
@test qc.measures_matrix == eye(4)

# Add measures
cr = ClassicalRegister(1)
setClassicalRegister!(qc, cr)
qc.measure([0], [0])
@test err(execute(backend, qc), [0.5, 0.5]) < 1e-5
@test qc.measures_matrix == [1 0 1 0; 0 1 0 1]

#############################################################################
n = 3
qr = QuantumRegister(7, "q")
cr = ClassicalRegister(3)
qc = QCircuit(qr, cr)
qr = [qr[1], qr[3], qr[5]]

qc.u3(qr)
for i in (n-2):-1:0
    i = i+1
    qc.cx(qr[i], qr[i+1])
    #qc.rzx(qr[i], qr[i+1])
end
qc.u3(qr)
for i in 0:(n-2)
    i = i+1
    qc.cx(qr[i], qr[i+1])
    #qc.rzx(qr[i], qr[i+1])
end
qc.u3(qr)
qc.measure([1, 3, 5], [0, 1, 2])

pa = [3.762539403251488, 0.001127890377849541, 2.1658488014045303, 6.174052504545824, 4.9629009414952305, 0.38255621799191375, 1.7274354859392205, 1.0857210633046896, 3.103638820041333, 3.758391851901853, 0.5765051731246599, 0.5798377641373272, 3.8759799027809825, 0.6188604600285638, 4.2218535852028385, 2.2390389868469462, 2.6674508252859623, 5.878739321230936, 1.340082725285039, 5.916066560675415, 3.625237411363865, 5.4397023221566725, 6.212355261857649, 2.5222234263059673, 5.8952667846227405, 6.113119466717698, 0.9810581564033299]

loss(params) = loss_expected_zero_state(execute(backend, qc, params))
dloss(params) = real(loss'(params))

ml = loss(pa)
mdl = dloss(pa)

sl, sdl = qderivative(qiskitBackend, qc, loss_expected_zero_state, pa, 4)
#println(join(pa, ", "))
@test abs(ml - sl) < 0.2
@test maximum(abs.(mdl - sdl)) < 2.0

# sl2, sdl2 = qderivative(qiskitBackend, ansact, loss_expected_zero_state, pa, 1)
# @test abs(ml - sl) - abs(ml - sl2) < 0.05
# @test maximum(abs.(mdl - sdl)) - maximum(abs.(mdl - sdl2)) < 0.1

ext: jl


text:
@doc doc"""
    MisclassLoss{R<:AbstractFloat} <: SupervisedLoss

Misclassification loss that assigns `1` for misclassified
examples and `0` otherwise. It is a generalization of
`ZeroOneLoss` for more than two classes.

The type parameter `R` specifies the result type of the
loss. Default type is double precision `R = Float64`.
"""
struct MisclassLoss{R<:AbstractFloat} <: SupervisedLoss end

MisclassLoss() = MisclassLoss{Float64}()

# specialize result type
result_type(loss::MisclassLoss{R}, t::Type, o::Type) where R = R

# return floating point to avoid big integers in aggregations
value(::MisclassLoss{R}, agreement::Bool) where R = ifelse(agreement, zero(R), one(R))
deriv(::MisclassLoss{R}, agreement::Bool) where R = zero(R)
deriv2(::MisclassLoss{R}, agreement::Bool) where R = zero(R)

value(loss::MisclassLoss, target, output) = value(loss, target == output)
deriv(loss::MisclassLoss, target, output) = deriv(loss, target == output)
deriv2(loss::MisclassLoss, target, output) = deriv2(loss, target == output)

# specialize aggregation methods (temporary workaround)
for FUN in (:value, :deriv, :deriv2)
    @eval begin
        # by default compute the element-wise result
        function ($FUN)(
                loss::MisclassLoss,
                targets::AbstractVector,
                outputs::AbstractVector)
            ($FUN)(loss, targets, outputs, AggMode.None())
        end

        function ($FUN)(
                loss::MisclassLoss,
                targets::AbstractVector,
                outputs::AbstractVector,
                ::AggMode.None)
            ($FUN).(loss, targets, outputs)
        end

        function ($FUN)(
                loss::MisclassLoss,
                targets::AbstractVector{Q},
                outputs::AbstractVector{T},
                ::AggMode.Sum) where {Q,T}
            out = zero(result_type(loss, Q, T))
            @inbounds @simd for I in CartesianIndices(size(outputs))
                out += ($FUN)(loss, targets[I], outputs[I])
            end
            out
        end

        function ($FUN)(
                loss::MisclassLoss,
                targets::AbstractVector{Q},
                outputs::AbstractVector{T},
                ::AggMode.Mean) where {Q,T}
            nrm = length(targets)
            out = zero(result_type(loss, Q, T))
            @inbounds @simd for I in CartesianIndices(size(outputs))
                out += ($FUN)(loss, targets[I], outputs[I])
            end
            out / nrm
        end

        function ($FUN)(
                loss::MisclassLoss,
                targets::AbstractVector{Q},
                outputs::AbstractVector{T},
                agg::AggMode.WeightedSum) where {Q,T}
            nrm = agg.normalize ? inv(sum(agg.weights)) : inv(one(sum(agg.weights)))
            out = zero(result_type(loss, Q, T)) * (agg.weights[1] * nrm)
            @inbounds @simd for I in CartesianIndices(size(outputs))
                out += ($FUN)(loss, targets[I], outputs[I]) * (agg.weights[I] * nrm)
            end
            out
        end

        function ($FUN)(
                loss::MisclassLoss,
                targets::AbstractVector{Q},
                outputs::AbstractVector{T},
                agg::AggMode.WeightedMean) where {Q,T}
            k   = length(outputs)
            nrm = agg.normalize ? inv(k * sum(agg.weights)) : inv(k * one(sum(agg.weights)))
            out = zero(result_type(loss, Q, T)) * (agg.weights[1] * nrm)
            @inbounds @simd for I in CartesianIndices(size(outputs))
                out += ($FUN)(loss, targets[I], outputs[I]) * (agg.weights[I] * nrm)
            end
            out
        end
    end
end

isminimizable(::MisclassLoss) = false
isdifferentiable(::MisclassLoss) = false
isdifferentiable(::MisclassLoss, at) = at != 0
istwicedifferentiable(::MisclassLoss) = false
istwicedifferentiable(::MisclassLoss, at) = at != 0
isnemitski(::MisclassLoss) = false
islipschitzcont(::MisclassLoss) = false
isconvex(::MisclassLoss) = false
isclasscalibrated(::MisclassLoss) = false
isclipable(::MisclassLoss) = false

# ===============================================================

@doc doc"""
    PoissonLoss <: SupervisedLoss

Loss under a Poisson noise distribution (KL-divergence)

``L(target, output) = exp(output) - target*output``
"""
struct PoissonLoss <: SupervisedLoss end

value(loss::PoissonLoss, target::Number, output::Number) = exp(output) - target*output
deriv(loss::PoissonLoss, target::Number, output::Number) = exp(output) - target
deriv2(loss::PoissonLoss, target::Number, output::Number) = exp(output)

isdifferentiable(::PoissonLoss) = true
isdifferentiable(::PoissonLoss, y, t) = true
istwicedifferentiable(::PoissonLoss) = true
istwicedifferentiable(::PoissonLoss, y, t) = true
islipschitzcont(::PoissonLoss) = false
isconvex(::PoissonLoss) = true
# TODO: isstrictlyconvex(::PoissonLoss) = ?
isstronglyconvex(::PoissonLoss) = false

# ===============================================================

@doc doc"""
    CrossEntropyLoss <: SupervisedLoss

Cross-entropy loss also known as log loss and logistic loss is defined as:

``L(target, output) = - target*ln(output) - (1-target)*ln(1-output)``
"""
struct CrossEntropyLoss <: SupervisedLoss end

function value(loss::CrossEntropyLoss, target::Number, output::Number)
    target >= 0 && target <=1 || error("target must be in [0,1]")
    output >= 0 && output <=1 || error("output must be in [0,1]")
    if target == 0
        -log(1 - output)
    elseif target == 1
        -log(output)
    else
        -(target * log(output) + (1-target) * log(1-output))
    end
end
deriv(loss::CrossEntropyLoss, target::Number, output::Number) = (1-target) / (1-output) - target / output
deriv2(loss::CrossEntropyLoss, target::Number, output::Number) = (1-target) / (1-output)^2 + target / output^2

isdifferentiable(::CrossEntropyLoss) = true
isdifferentiable(::CrossEntropyLoss, y, t) = t != 0 && t != 1
istwicedifferentiable(::CrossEntropyLoss) = true
istwicedifferentiable(::CrossEntropyLoss, y, t) = t != 0 && t != 1
isconvex(::CrossEntropyLoss) = true

# ===============================================================
# L(target, output) = sign(agreement) < 0 ? 1 : 0


text:
using NCDatasets

const SOCAT_FILE = "SOCATv5_all.tsv"::String
const DECORRELATION_LENGTHS_FILE = "decorrelation_lengths.csv"::String
const CELL_SIZE = 2.5::Float64

# Get the longitudes
function getlongitudes()
    local longitudes::Array{Float64, 1} = Array{Float64, 1}(convert(Int64, 360 / CELL_SIZE))
    for i in 1:length(longitudes)
        longitudes[i] = i * CELL_SIZE - (CELL_SIZE / 2)
    end

    return longitudes
end

# Get the latitudes
function getlatitudes()
    local latitudes::Array{Float64, 1} = Array{Float64, 1}(convert(Int64, 180 / CELL_SIZE))
    for i in 1:length(latitudes)
        latitudes[i] = i * CELL_SIZE - (CELL_SIZE / 2) - 90
    end

    return latitudes
end

# Calculate the cell indices from a longitude and latitude
function getcellindex(longitude::Float64, latitude::Float64)
    local loncell::Int64 = floor(longitude / CELL_SIZE) + 1
    local latcell::Int64 = floor((latitude + 90) / CELL_SIZE) + 1

    if loncell == 145
        loncell = 1
    end

    return loncell, latcell
end

# Add the specified decorrelation length to cells specified
function assigndatasetcells(totals::Array{Int64, 2}, counts::Array{Int64, 2}, cellsindataset::Array{Bool, 2}, decorrelationlength::Int64)
    if decorrelationlength != -1
        for lon in 1:size(cellsindataset)[1]
            for lat in 1:size(cellsindataset)[2]

                if cellsindataset[lon, lat]
                    totals[lon, lat] = totals[lon, lat] + decorrelationlength
                    counts[lon, lat] = counts[lon, lat] + 1
                end
            end
        end
    end
end

# Main function
function run()

    # Initialise data structures
    local totals::Array{Int64, 2} = zeros(convert(Int64, 360 / CELL_SIZE), convert(Int64, 180 / CELL_SIZE))
    local counts::Array{Int64, 2} = zeros(convert(Int64, 360 / CELL_SIZE), convert(Int64, 180 / CELL_SIZE))

    # Read in the decorrelation lengths
    local decorrelationlengths::Array{String, 2} = readdlm(DECORRELATION_LENGTHS_FILE, ',', String, '\n', header=false)

    local currentdataset::String = ""
    local currentdecorrelationlength::Int64 = -1
    local datasetcells::Array{Bool, 2} = falses(convert(Int64, 360 / CELL_SIZE), convert(Int64, 180 / CELL_SIZE))

    local socatchan::IOStream = open(SOCAT_FILE)

    currentline = readline(socatchan)
    local linecount = 1

    # Loop until the end of the file
    while length(currentline) > 0

        fields = split(currentline, "\t")
        dataset = fields[1]

        if (dataset != currentdataset)
            if !isnull(currentdataset)
                assigndatasetcells(totals, counts, datasetcells, currentdecorrelationlength)
            end

            currentdataset = dataset
            rowsearch = find(decorrelationlengths[:,1] .== "$(dataset).tsv")
            if size(rowsearch) == 0 || size(rowsearch)[1] == 0
                decorrelationrow = -1
            else
                decorrelationrow = rowsearch[1]
            end

            if decorrelationrow != -1
                currentdecorrelationlength = parse(Int, decorrelationlengths[decorrelationrow[1], 2])
            else
                # The data set doesn't have a decorrelation length,
                # so we'll ignore it later on
                currentdecorrelationlength = -1
            end

            datasetcells = falses(convert(Int64, 360 / CELL_SIZE), convert(Int64, 180 / CELL_SIZE))

            print("\033[1K\r$dataset ($linecount)")
        end

        lon = parse(Float64, fields[11])
        lat = parse(Float64, fields[12])

        cell = getcellindex(lon, lat)
        datasetcells[cell[1], cell[2]] = true

        currentline = readline(socatchan)
        linecount += 1
    end

    if !isnull(currentdataset)
        assigndatasetcells(totals, counts, datasetcells, currentdecorrelationlength)
    end

    # Calculate the mean decorrelation lengths and write them to disk
    print("\033[1K\rWriting to disk...")
    local meandecorrelationlengths::Array{Float64, 2} = fill(NaN, (convert(Int64, 360 / CELL_SIZE), convert(Int64, 180 / CELL_SIZE)))
    for i in 1:convert(Int64, 360 / CELL_SIZE)
        for j in 1:convert(Int64, 180 / CELL_SIZE)
            if counts[i, j] > 0
                meandecorrelationlengths[i, j] = totals[i, j] / counts[i, j]
            end
        end
    end

    ds::Dataset = Dataset("mean_decorrelation_lengths.nc", "c")

    defDim(ds, "longitude", convert(Int64, 360 / CELL_SIZE))
    defDim(ds, "latitude", convert(Int64, 180 / CELL_SIZE))

    local lonvar::NCDatasets.CFVariable = defVar(ds, "longitude", Float32, ["longitude"])
    lonvar.attrib["units"] = "degrees_east"
    local latvar::NCDatasets.CFVariable = defVar(ds, "latitude", Float32, ["latitude"])
    latvar.attrib["units"] = "degrees_north"

    local meanvar::NCDatasets.CFVariable = defVar(ds, "decorrelation_lengths", Float64, ["longitude", "latitude"], fillvalue=-999.9)
    local totalvar::NCDatasets.CFVariable = defVar(ds, "totals", Float64, ["longitude", "latitude"], fillvalue=-999.9)
    local countvar::NCDatasets.CFVariable = defVar(ds, "counts", Float64, ["longitude", "latitude"], fillvalue=-999.9)

    lonvar[:] = getlongitudes()
    latvar[:] = getlatitudes()
    meanvar[:, :] = meandecorrelationlengths
    totalvar[:, :] = totals
    countvar[:, :] = counts

    close(ds)

    println()
end

run()


text:
export AdiabaticIdealGasEOS, MONATOMIC, DIATOMIC, PHYSICAL, NORMALIZED

abstract type Atomicity end
struct Monatomic <: Atomicity end
struct Diatomic <: Atomicity end
const MONATOMIC = Monatomic()
const DIATOMIC = Diatomic()

abstract type EOSType end
struct Physical <: EOSType end
struct Normalized <: EOSType end
const PHYSICAL = Physical()
const NORMALIZED = Normalized()

degrees_of_freedom(::Dim{N}, ::Monatomic) where {N} = N
degrees_of_freedom(::Dim{N}, ::Diatomic) where {N} = N + 2

"Equation of state for an ideal adiabatic gas."
struct AdiabaticIdealGasEOS <: EquationOfState
    "Adiabatic index."
    γ::Float
    "Pressure from internal energy density."
    P_from_e::Float
    "Internal energy density from mass density and temperature."
    e_from_ρT::Float

    """
        AdiabaticIdealGasEOS(
            μ::Number,
            ::Physical;
            dim::Dim = DIM3,
            atomicity::Atomicity = MONATOMIC,
            units::UnitSystem = SI,
        )

    Create a new ideal equation of state for an adiabatic gas with the given
    mean molecular mass `μ`, dimensionality, atomicity and unit system.
    """
    function AdiabaticIdealGasEOS(
        μ::Number,
        ::Physical;
        dim::Dim = DIM3,
        atomicity::Atomicity = MONATOMIC,
        units::UnitSystem = SI,
    )
        @check μ > 0
        γ = 1 + 2 / degrees_of_freedom(dim, atomicity)
        P_from_e = γ - 1
        e_from_ρT = units.k / (μ * P_from_e)
        new(γ, P_from_e, e_from_ρT)
    end

    """
        AdiabaticIdealGasEOS(γ::Number, ::Normalized)

    Create a new ideal equation of state, normalized so `P = ρT`, for an
    adiabatic gas with the given adiabatic index `γ`.
    """
    function AdiabaticIdealGasEOS(γ::Number, ::Normalized)
        @check γ > 0
        P_from_e = γ - 1
        e_from_ρT = 1 / P_from_e
        new(γ, P_from_e, e_from_ρT)
    end
end

energy_density_from(eos::AdiabaticIdealGasEOS, ρ::Number, u::Number) =
    @dbgassert(ρ ≥ 0) && @dbgassert(u ≥ 0) && ρ * u

gas_pressure_from(eos::AdiabaticIdealGasEOS, ρ::Number, u::Number) =
    eos.P_from_e * energy_density_from(eos, ρ, u)

temperature_from(eos::AdiabaticIdealGasEOS, u::Number) =
    @dbgassert(u ≥ 0) && u / eos.e_from_ρT

mass_density_from(eos::AdiabaticIdealGasEOS, P::Number, u::Number) =
    @dbgassert(P ≥ 0) && @dbgassert(u > 0) && P / (eos.P_from_e * u)

specific_energy_from(eos::AdiabaticIdealGasEOS, T::Number) =
    @dbgassert(T ≥ 0) && eos.e_from_ρT * T

specific_energy_from(eos::AdiabaticIdealGasEOS, ρ::Number, P::Number) =
    @dbgassert(P ≥ 0) && @dbgassert(ρ > 0) && P / (eos.P_from_e * ρ)

soundspeed²_from(eos::AdiabaticIdealGasEOS, u::Number) =
    @dbgassert(u ≥ 0) && eos.γ * eos.P_from_e * u

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# compiler support for working with run-time libraries

function link_library!(ctx::CompilerContext, mod::LLVM.Module, lib::LLVM.Module)
    # linking is destructive, so copy the library
    lib = LLVM.Module(lib)

    # save list of external functions
    exports = String[]
    for f in functions(mod)
        fn = LLVM.name(f)
        if !haskey(functions(lib), fn)
            push!(exports, fn)
        end
    end

    link!(mod, lib)

    ModulePassManager() do pm
        # internalize all functions that aren't exports
        internalize!(pm, exports)

        # eliminate all unused internal functions
        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)

        run!(pm, mod)
    end
end

const libcache = Dict{String, LLVM.Module}()

#
# CUDA device library
#

function find_libdevice(cap)
    CUDAnative.configured || return
    global libdevice

    if isa(libdevice, Dict)
        # select the most recent & compatible library
        vers = keys(CUDAnative.libdevice)
        compat_vers = Set(ver for ver in vers if ver <= cap)
        isempty(compat_vers) && error("No compatible CUDA device library available")
        ver = maximum(compat_vers)
        path = libdevice[ver]
    else
        libdevice
    end
end

function load_libdevice(cap)
    path = find_libdevice(cap)

    get!(libcache, path) do
        open(path) do io
            parse(LLVM.Module, read(path), JuliaContext())
        end
    end
end

function link_libdevice!(ctx::CompilerContext, mod::LLVM.Module, lib::LLVM.Module)
    # override libdevice's triple and datalayout to avoid warnings
    triple!(lib, triple(mod))
    datalayout!(lib, datalayout(mod))

    link_library!(ctx, mod, lib)

    ModulePassManager() do pm
        push!(metadata(mod), "nvvm-reflect-ftz",
              MDNode([ConstantInt(Int32(1), JuliaContext())]))
        # TODO: run the reflect pass?
        run!(pm, mod)
    end
end

#
# CUDAnative run-time library
#

# remove existing runtime libraries globally,
# so any change to CUDAnative triggers recompilation
rm(joinpath(@__DIR__, "..", "..", "deps", "runtime"); recursive=true, force=true)

## higher-level functionality to work with runtime functions

function LLVM.call!(builder, rt::Runtime.RuntimeMethodInstance, args=LLVM.Value[])
    bb = position(builder)
    f = LLVM.parent(bb)
    mod = LLVM.parent(f)

    # get or create a function prototype
    if haskey(functions(mod), rt.llvm_name)
        f = functions(mod)[rt.llvm_name]
        ft = eltype(llvmtype(f))
    else
        ft = LLVM.FunctionType(rt.llvm_return_type, rt.llvm_types)
        f = LLVM.Function(mod, rt.llvm_name, ft)
    end

    # runtime functions are written in Julia, while we're calling from LLVM,
    # this often results in argument type mismatches. try to fix some here.
    for (i,arg) in enumerate(args)
        if llvmtype(arg) != parameters(ft)[i]
            if (llvmtype(arg) isa LLVM.PointerType) &&
               (parameters(ft)[i] isa LLVM.IntegerType)
                # Julia pointers are passed as integers
                args[i] = ptrtoint!(builder, args[i], parameters(ft)[i])
            else
                error("Don't know how to convert ", arg, " argument to ", parameters(ft)[i])
            end
        end
    end

    call!(builder, f, args)
end

## functionality to build the runtime library

function emit_function!(mod, cap, f, types, name)
    tt = Base.to_tuple_type(types)
    ctx = CompilerContext(f, tt, cap, #= kernel =# false)
    new_mod, entry = irgen(ctx)
    entry = optimize!(ctx, new_mod, entry)
    LLVM.name!(entry, name)

    link!(mod, new_mod)
end

function build_runtime(cap)
    mod = LLVM.Module("CUDAnative run-time library", JuliaContext())

    for method in values(Runtime.methods)
        emit_function!(mod, cap, method.def, method.types, method.llvm_name)
    end

    mod
end

function load_runtime(cap)
    name = "cudanative.$(cap.major)$(cap.minor).bc"
    path = joinpath(@__DIR__, "..", "..", "deps", "runtime", name)
    mkpath(dirname(path))

    get!(libcache, path) do
        if ispath(path)
            open(path) do io
                parse(LLVM.Module, read(io), JuliaContext())
            end
        else
            @info "Building the CUDAnative run-time library for your sm_$(cap.major)$(cap.minor) device, this might take a while..."
            lib = build_runtime(cap)
            open(path, "w") do io
                write(io, lib)
            end
            lib
        end
    end
end

max_stars_repo_name: Cyvadra/CUDAnative

max_issues_repo_name: Cyvadra/CUDAnative

max_forks_repo_name: Cyvadra/CUDAnative


text:
# small scrip to compute the solution of Lippmann-Schwinger equation
# We test the types introduced in FastConvolution.jl
# we test that the application is fast and that the construction
# is performed fast.

using PyPlot
using IterativeSolvers

include("../src/FastConvolution.jl")
include("../src/Preconditioner.jl")

# setting the number of threads for the FFT and BLAS
# libraries (please set them to match the number of
# physical cores in your system)
FFTW.set_num_threads(4);
BLAS.set_num_threads(4);

#Defining Omega
#Defining Omega
h = 1/32
k = 2.0

# size of box
a  = 1
x = collect(-a/2:h:a/2-h)
y = collect(-a/2:h:a/2-h)
z = collect(-a/2:h:a/2-h)

(n,m,l) = length(x), length(y), length(z)

N = n*m*l

# creating the 3D tensor grid
# TODO: write a
X = zeros(n,m,l); Y = zeros(n,m,l); Z = zeros(n,m,l);
for i=1:n, j=1:m, p=1:l
    X[i,j,p] =  x[i];
    Y[i,j,p] =  y[j];
    Z[i,j,p] =  z[p];
end

X = X[:]; Y = Y[:]; Z = Z[:];

# Defining the smooth perturbation of the slowness
sigma = 0.05
renorm_factor = 1/(sigma^3*(2*pi)^(3/2));
rho(x,y,z) = renorm_factor*exp((-1/(2*sigma^2))*(x.^2 + y.^2 + z.^2));

Rho = rho(X,Y,Z) + 0*1im;
RHO = reshape(Rho, n,m,l);

figure(5); clf();
imshow(real(RHO[:,:,round(Integer, end/2)]))
## You can choose between Duan Rohklin trapezoidal quadrature
# fastconv = buildFastConvolution(x,y,h,k,nu)

# defining a dummy nu
nu(x,y,z) = 1*exp(-40*(x.^2 + y.^2 + z.^2)).*(abs(x).<0.48).*(abs(y).<0.48).*(abs(z).<0.48);

# or Greengard Vico Quadrature (this is not optimized and is 2-3 times slower)
fastconv = buildFastConvolution3D(x,y,z,X,Y,Z,h,k,nu, quadRule = "Greengard_Vico");

u_rho = FFTconvolution(fastconv, Rho);

solRef(x,y,z) = exp(-sigma^2*k^2/2)./(4*pi*sqrt(x.^2+y.^2+z.^2)).*
                (real(exp(-1im*k*sqrt(x.^2+y.^2+z.^2)).*erf((2*sigma^2*1im*k - 2*sqrt(x.^2 + y.^2 + z.^2))/(2*sqrt(2*sigma^2)))) - 1im*sin(k*sqrt(x.^2 + y.^2 + z.^2)))

uRef = solRef(X,Y,Z);

# plotting the solution
figure(1)
clf()
UU = reshape(u_rho,n,m,l)
imshow(real(UU[:,:,round(Integer, end/2)])); colorbar();

figure(2)
clf();
UURef = reshape(uRef ,n,m,l)
imshow(real(UURef[:,:,round(Integer, end/2)])); colorbar();

figure(3)
clf();
imshow(abs((UURef-UU)[:,:,round(Integer, end/2)])); colorbar()
println("Error in L infinity norm is ", maximum(abs(u_rho - uRef) ))

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
###################################################################################################
# FUNCTIONS
###################################################################################################
"""
    `generate_xcal(N)`

Generate state space 𝒳 for the methylation vector with N CpG sites.

# Examples
```julia-repl
julia> xcal=CpelAsm.generate_xcal(2)
4-element Array{Vector{Int64},1}:
 [-1, -1]
 [1, -1]
 [-1, 1]
 [1, 1]
```
"""
function generate_xcal(n::Int64)::Array{Vector{Int64},1}

    # Generate iterative object
    xcal = Vector{Int64}[]
    [push!(xcal,2 .* digits(i, base=2, pad=n) - ones(n)) for i in 0:(2^n-1)]

    # Return 𝒳
    return xcal

end # end generate_xcal
"""
    `gen_ising_full_data(M,N;α,β)`

Generate M FULLY observed reads with [N1,...,NK] CpG sites from energy function
U(X) parametrized by [α1,...,αK] and β.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=CpelAsm.gen_ising_full_data(10,[4])
10-element Array{Array{Int64,1},1}:
 [-1, 1, -1, -1]
 [1, -1, 1, -1]
 [-1, 1, 1, -1]
 [1, 1, 1, -1]
 [1, 1, 1, -1]
 [-1, 1, -1, 1]
 [-1, 1, -1, 1]
 [-1, 1, -1, 1]
 [-1, -1, 1, 1]
 [1, 1, 1, 1]
```
"""
function gen_ising_full_data(m::Int64,n::Vector{Int64};a=zeros(length(n)),b=0.0)::Array{Vector{Int64},1}

    # Generate iterative object
    xcal = generate_xcal(sum(n))

    # Create function for energy computation
    Ux_fun = create_Ux(n,a,b)

    # Compute probability over state space
    p = zeros(2^sum(n))
    p = [exp(-Ux_fun(xcal[k])) for k in 1:length(xcal)]
    p = p ./ sum(p)

    # Sample full observations from model
    indices = rand(Distributions.Multinomial(m, p))
    indices = vcat([repeat([i],indices[i]) for i in 1:length(indices)]...)

    # Return observations
    return xcal[indices]

end # end gen_ising_full_data
"""
    `gen_ising_part_data(M,N;α,β)`

Generate M partially observed reads with [N1,...,NK] CpG sites from energy function U(X)
parametrized by [α1,...,αK] and β. The values code for:
    1:  Methylated CpG site.
    -1: Unmethylated CpG site.
    0:  Unobserved CpG site.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=gen_ising_part_data(10,[4])
10-element Array{Vector{Int64},1}:
[-1, -1, -1, -1]
[1, -1, -1, -1]
[-1, 1, -1, -1]
[-1, -1, 1, -1]
[1, -1, 1, -1]
[1, -1, -1, 1]
[0, 1, -1, 1]
[1, -1, 1, 1]
[1, -1, 1, 1]
[-1, 1, 1, 1]
```
"""
function gen_ising_part_data(m::Int64,n::Vector{Int64};a=zeros(length(n)),b=0.0)::Array{Vector{Int64},1}

    # Get full data
    full_data = gen_ising_full_data(m,n;a=a,b=b)

    # Add missing values on sides
    missleft = rand(Distributions.Binomial(sum(n), 0.01),m)
    missright = rand(Distributions.Binomial(sum(n), 0.01),m)
    for i in 1:m
        full_data[i][1:missleft[i]] = fill(0, missleft[i])
        full_data[i][(sum(n)-missright[i]+1):end] = fill(0, missright[i])
    end

    # Return partial observations
    return full_data

end # end gen_ising_part_data
"""
    `gen_mult_full_data(M;N,p)`

Generate M FULLY observed reads from a multinomial model with vector of probabilities
`p=[p_1,...,p_{2^N}]`.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=CpelAsm.gen_mult_full_data(2)
2-element Array{Array{Int64,1},1}:
 [1, 1, -1, 1]
 [1, 1, -1, 1]
```
"""
function gen_mult_full_data(m::Int64;n=4,p=1.0/(2.0^n)*ones(2^n))::Array{Vector{Int64},1}

    # Generate iterative object
    xcal = generate_xcal(n)

    # Sample full observations from model
    indices = rand(Distributions.Multinomial(m, p))
    indices = vcat([repeat([i],indices[i]) for i in 1:length(indices)]...)

    # Return observations
    return xcal[indices]

end # end gen_mult_full_data
"""
    `gen_mult_part_data(M;N,p)`

Generate M partially observed reads with N CpG sites from a multinomial model with vector
probabilities `p=[p_1,...,p_{2^N}]`. The values code for:
    1:  Methylated CpG site.
    -1: Unmethylated CpG site.
    0:  Unobserved CpG site.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=CpelAsm.gen_mult_part_data(10)
10-element Array{Array{Int64,1},1}:
 [0, 0, 0, 0]
 [1, -1, 1, -1]
 [0, 1, 1, -1]
 [1, 1, 1, 0]
 [1, 1, 1, 0]
 [0, 1, -1, 1]
 [0, 1, -1, 1]
 [0, 1, -1, 1]
 [0, -1, 1, 1]
 [1, 1, 1, 1]
```
"""
function gen_mult_part_data(m::Int64;n=4,p=1.0/(2.0^n)*ones(2^n))::Array{Vector{Int64},1}

    # Get full data
    full_data = gen_mult_full_data(m;n=n,p=p)

    # Add missing values on sides
    missleft = rand(Distributions.Binomial(n+1, 0.05),m)
    missright = rand(Distributions.Binomial(n+1, 0.05),m)
    for i in 1:m
        full_data[i][1:missleft[i]] = fill(0, missleft[i])
        full_data[i][(n-missright[i]+1):end] = fill(0, missright[i])
    end

    # Return partial observations
    return full_data

end # end gen_mult_part_data


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.3773826343471955 <= q <= 0.8255956851677421)

                     add_NL_constraint(m, :(swish(-0.5817410001457377 + -0.38051045291676777*swish(-0.7752766980375418 + 0.20887750352529233*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.828517562717189*swish(-0.31567568822764036 + 0.26872838410166544*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.13996198025628548*swish(0.6085123540494801 + -0.012952465394535206*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) + swish(-0.08252400134510385 + 0.3442818429504264*swish(-0.7752766980375418 + 0.20887750352529233*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.7813134096358496*swish(-0.31567568822764036 + 0.26872838410166544*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.5888247957791601*swish(0.6085123540494801 + -0.012952465394535206*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) + swish(-0.2140493400207517 + 0.6401317985056014*swish(-0.7752766980375418 + 0.20887750352529233*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.4453350709432602*swish(-0.31567568822764036 + 0.26872838410166544*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + 0.9148591417614278*swish(0.6085123540494801 + -0.012952465394535206*swish(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*swish(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*swish(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function compute_ν(n, dim)
    if dim == 2
        ν = n^2
    else
        ν = n * (n + 1)
    end

    ν
end


text:
module Move_FourBar_noMacros

import Modia3D
using  Modia3D.ModiaMath
using  Modia3D.Unitful

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
massProperties = "Aluminium"
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

# Base dimensions
length1  = 1.0
width    = 0.2
height   = 0.2
diameter = 0.1
groundWidth  = 1.0
groundHeight = 0.01
cyl  = Modia3D.Cylinder(width/2, 1.2*width; material=vmat2)
Tcyl = ModiaMath.rot1(90u"°")
TGround = ModiaMath.rot1(-90u"°")

# -----------------------------------------------------------------------
# @assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=Ly) begin
#   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vmat1))
#   frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
#   frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])
#   cyl1   = Modia3D.Object3D(frame1, cyl)
#   cyl2   = Modia3D.Object3D(frame2, cyl)
# end
mutable struct Bar <: Modia3D.AbstractAssembly
   _internal::Modia3D.AssemblyInternal
   Lx
   Ly
   Lz
   frame0
   frame1
   frame2
   cyl1
   cyl2

   function Bar(;sceneOptions=nothing,Lx = 0.1, Ly=Lx/5, Lz=Ly)
      this = new( Modia3D.AssemblyInternal(:Bar,sceneOptions) )
      this.Lx = Lx
      this.Ly = Ly
      this.Lz = Lz

      frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vmat1))
      Modia3D.initAssemblyComponent!(this, frame0, "frame0")

      frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
      Modia3D.initAssemblyComponent!(this, frame1, "frame1")

      frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])
      Modia3D.initAssemblyComponent!(this, frame2, "frame2")

      cyl1   = Modia3D.Object3D(frame1, cyl)
      Modia3D.initAssemblyComponent!(this, cyl1, "cyl1")

      cyl2   = Modia3D.Object3D(frame2, cyl)
      Modia3D.initAssemblyComponent!(this, cyl2, "cyl2")

      return this
   end
end
# -----------------------------------------------------------------------

# -----------------------------------------------------------------------
# @assembly Fourbar begin
#    world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
#    frame1 = Modia3D.Object3D(world ; r=[0.5*length1, 0.0, groundWidth/2])
#    frame2 = Modia3D.Object3D(frame1; r=[length1    , 0.0, 0.0])
#    ground = Modia3D.Object3D(world , Modia3D.Box(3*length1, groundWidth, groundHeight; material=groundMaterial);
#                                               r = [1.5*length1, -groundHeight/2, groundWidth/2],
#                                               R = Modia3D.rot(;n1 = [-1,0,0], n2 = [0,0,1]))
#    bar1   = Bar(Lx=length1)
#    bar2   = Bar(Lx=length1)
#    bar3   = Bar(Lx=length1)
#
#    rev1   = Modia3D.Revolute(frame1     , bar1.frame1; phi =  pi/2)
#    rev2   = Modia3D.Revolute(bar1.frame2, bar2.frame1; phi = -pi/2)
#    rev3   = Modia3D.Revolute(frame2     , bar3.frame1; phi =  pi/2)
#    rev4   = Modia3D.Revolute(bar3.frame2, bar2.frame2; phi = -pi/2)
# end
mutable struct Fourbar <: Modia3D.AbstractAssembly
   _internal::Modia3D.AssemblyInternal
   world
   frame1
   frame2
   ground
   bar1
   bar2
   bar3
   rev1
   rev2
   rev3
   rev4

   function Fourbar(;sceneOptions=nothing)
      this = new( Modia3D.AssemblyInternal(:Fourbar, sceneOptions) )

      world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
      Modia3D.initAssemblyComponent!(this, world , "world")

      frame1 = Modia3D.Object3D(world ; r=[0.5*length1, 0.0, groundWidth/2])
      Modia3D.initAssemblyComponent!(this, frame1, "frame1")

      frame2 = Modia3D.Object3D(frame1; r=[length1, 0.0, 0.0])
      Modia3D.initAssemblyComponent!(this, frame2, "frame2")

      ground = Modia3D.Object3D(world , Modia3D.Box(3*length1, groundWidth, groundHeight; material=groundMaterial);
                                                 r = [1.5*length1, -groundHeight/2, groundWidth/2],
                                                 R = ModiaMath.rot_nxy([-1,0,0], [0,0,1]))
      Modia3D.initAssemblyComponent!(this, ground, "ground")

      bar1 = Bar(Lx=length1)
      Modia3D.initAssemblyComponent!(this, bar1, "bar1")

      bar2 = Bar(Lx=length1)
      Modia3D.initAssemblyComponent!(this, bar2, "bar2")

      bar3 = Bar(Lx=length1)
      Modia3D.initAssemblyComponent!(this, bar3, "bar3")

      rev1 = Modia3D.Revolute(frame1, bar1.frame1; phi_start =  pi/2)
      Modia3D.initAssemblyComponent!(this, rev1, "rev1")

      rev2 = Modia3D.Revolute(bar1.frame2, bar2.frame1; phi_start = -pi/2)
      Modia3D.initAssemblyComponent!(this, rev2, "rev2")

      rev3 = Modia3D.Revolute(frame2, bar3.frame1; phi_start =  pi/2)
      Modia3D.initAssemblyComponent!(this, rev3, "rev3")

      rev4 = Modia3D.Revolute(bar3.frame2, bar2.frame2; phi_start = -pi/2)
      Modia3D.initAssemblyComponent!(this, rev4, "rev4")

      return this
   end
end
# -----------------------------------------------------------------------

# Kinematic simulation
fourbar = Fourbar(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,
                                                    defaultFrameLength=0.3,
                                                    enableContactDetection=false))
Modia3D.initAnalysis!(fourbar)
tStart=0.0
tEnd  =1.0
Modia3D.driveJoint!(fourbar.rev1)  # Define that rev1 is driven (setAngle! can be called on it)
Modia3D.driveJoint!(fourbar.rev2)  # Define that rev2 is driven (setAngle! can be called on it)

for time = range(tStart, stop=tEnd, length=101)
  # update positional degrees of freedom
  delta_phi = Modia3D.linearMovement(pi/3, tStart, tEnd, time)
  Modia3D.setAngle!(fourbar.rev1,  pi/2 - delta_phi)
  Modia3D.setAngle!(fourbar.rev2, -pi/2 + delta_phi)

  Modia3D.updatePosition!(fourbar)

  # Visualize
  Modia3D.visualize!(fourbar,time)
end

Modia3D.closeAnalysis!(fourbar)

println("... success of Move_FourBar_noMacros.jl!")
end


text:
# Julia typed form of the Harrison solid state table
# Walter A. Harrison - Electronic structure & the property of solids

export SolidStateTableEntry, SSTE, SolidStateTable, SST, ShiSST
export HarrisonIME, HarrisonShiAltIME, ShiIME

# Interatomic matrix elements (eta, η)
# Harrison, 'Solid State Table' parameters, top right table
# TODO: Check; and wrap these in a struct
struct InterAtomic
    ηssσ::Float64
    ηspσ::Float64
    ηppσ::Float64
    ηppπ::Float64

    ηsdσ::Float64
    ηpdσ::Float64
    ηpdπ::Float64

    ηddσ::Float64
    ηddπ::Float64
    ηddδ::Float64
end

# Nb: These are from the Dover reprint of Harrison's book (El-strucut & the prop of solids)
HarrisonIME=InterAtomic(-1.40,+1.84,+3.24,-0.81, -3.16,-2.95,+1.36, -16.2,+8.75,+0.0)

# Shi and Papaconstantopoulos give quite different values for Harrison;
# s and p are totally different; d mostly the same except for ddδ (last value)
#   Table I https://doi.org/ 10.1103/PhysRevB.70.205101
HarrisonShiAltIME=InterAtomic(-1.32,+1.42,+2.22,-0.63, -3.16,-2.95,+1.36, -16.2,+8.75,-2.39)

# Shi and Papconstantopoulos - improved values for d-block
#   See Table I https://doi.org/ 10.1103/PhysRevB.70.205101
ShiIME=     InterAtomic(-0.90,+1.44,+2.19,-0.03, -3.12,-4.26,+2.08, -21.22,+12.60,-2.29)

# Vaguely structured to look like the Solid State Table of the Elements
# Using struct so type immutable; allows for efficient memory use and type stability
# Current just considering 'simple atoms' form - SP hybridising semiconductors
struct SolidStateTableEntry
    name::AbstractString
    Z::Int

    ϵs::Float64
    ϵp::Float64
    ϵd::Float64

    kF::Float64
    rc::Float64 # also re-used for r0 (or ro?) of d-block
    ri::Float64

    γ::Float64

    mass::Float64
end

SSTE=SolidStateTableEntry # an alias to shorten the following lines; the full name appears in the REPL etc.

# Constructor methods
# SIMPLE ATOMS
SSTp(name,Z,ϵs,ϵp,kF,rc,ri,mass)=SSTE(name,Z,ϵs,ϵp,NaN, kF,rc,ri, NaN, mass)
# TRANSITION METALS (we call this d for 'd block')
SSTd(name,Z,ϵd,rd,kd,ro,ri,mass)=SSTE(name,Z,NaN,NaN,ϵd, NaN, ro, ri, NaN, mass)
# F-SHELL METALS
SSTf()=SSTE() # TODO: Implement f-shell in data structure

# Shi-Papaconstantopoulos
# 2004 https://doi.org/10.1103/PhysRevB.70.205101 ; 21st century improved d-orbital values
SSTShi(name,Z,ϵs,ϵp,ϵd,γ,rd,mass)=SSTE(name,Z, ϵs,ϵp,ϵd, NaN, NaN, rd, γ, mass)

# OK; here is the actual data - as in Harrison's book
SolidStateTable=Array{SSTE}(92) # only up to Z=92 , uranium currently!
SST=SolidStateTable

ShiSST=Array{SSTE}(92)
ShiSST[22]=SSTShi("Ti",22, 0.51252, 0.79759, 0.21879, 0.92307, 1.85267, 47.90)

#  "If it's not checked, it's wrong" - Slogan at Bletchley Park
# These values are not well checked.
SST[1]=SSTp("H" ,1,   NaN,NaN,      NaN,NaN,NaN,      1.0)
SST[2]=SSTp("He",2,   23.4,NaN,     NaN,NaN,NaN,      4.0)
SST[3]=SSTp("Li",3,   5.48,NaN,     1.13,0.92,0.68,   6.94)

SST[4]=SSTp("Be",4,   8.17,4.14,    1.94,0.58,0.30,   9.01)
SST[5]=SSTp("B", 5,   12.54,6.64,   NaN,0.44,0.16,    10.81)
SST[6]=SSTp("C", 6,   17.52,8.97,   2.76,0.37,NaN,    12.01)
SST[7]=SSTp("N", 7,   23.04,11.47,  NaN,NaN,NaN,      14.01)
SST[8]=SSTp("O", 8,   29.14,14.13,  0,0.42,1.46,      16.00)
# OK, not trying to be comprehensive at this point - let's just get some tetrahedral semiconductors on!

SST[12]=SSTp("Mg",12, 6.86, 2.99,   1.37,0.74,0.65,   24.31)

SST[14]=SSTp("Si",14, 13.55,6.52,   1.81,0.56,0.38,   28.09)
SST[15]=SSTp("P", 15, 17.10,8.33,   NaN,0.51,NaN,     30.97)
SST[16]=SSTp("S", 16, 20.80,10.27,  NaN,0.47,1.0,     32.06)

SST[22]=SSTd("Ti",22, 11.04, 1.08,  1.17, 1.61,0.90,   47.90) # TODO: Check d-block.

# Cu is special - also has d-block params?
SST[29]=SSTp("Cu",29, 6.92,1.83,    1.36,NaN,NaN,     63.54) # 'simple' atom form
SST[30]=SSTp("Zn",30, 8.40,3.38,    1.59,0.59,NaN,    65.37)
SST[31]=SSTp("Ga",31, 11.37,4.90,   1.66,0.59,NaN,    69.82)
SST[32]=SSTp("Ge",32, 14.38,6.36,   1.74,0.54,NaN,    72.59)
SST[33]=SSTp("As",33, 17.33,7.91,   NaN,0.51, NaN,    74.92)
SST[34]=SSTp("Se",34, 20.32,9.53,   NaN,0.50, NaN,    78.96)

SST[50]=SSTp("Sn",50, 12.50,5.94,   1.63,0.59,NaN,    118.7)

SST[53]=SSTp("I" ,53, 19.42,9.97,   NaN,NaN,2.16,     126.9)
SST[82]=SSTp("Pb",82, 12.07,5.77,   1.58,0.57,NaN,    207.2)