text:
Gm, Gn = 46, 46;
Gcolptr = [1, 4, 7, 11, 15, 18, 21, 25, 30, 34, 37, 41, 46, 50, 56, 60, 64, 69, 73, 77, 81, 86, 90, 93, 98, 102, 106, 110, 114, 117, 121, 126, 130, 134, 139, 143, 147, 152, 157, 161, 165, 169, 172, 175, 179, 186, 189];
Growval = [1, 2, 6, 1, 2, 3, 2, 3, 4, 7, 3, 4, 5, 8, 4, 5, 9, 1, 6, 10, 3, 7, 8, 11, 4, 7, 8, 9, 12, 5, 8, 9, 13, 6, 10, 14, 7, 11, 12, 16, 8, 11, 12, 13, 17, 9, 12, 13, 18, 10, 14, 15, 19, 24, 30, 14, 15, 16, 19, 11, 15, 16, 17, 12, 16, 17, 18, 20, 13, 17, 18, 21, 14, 15, 19, 24, 17, 20, 21, 26, 18, 20, 21, 22, 27, 21, 22, 23, 28, 22, 23, 29, 14, 19, 24, 25, 31, 24, 25, 26, 32, 20, 25, 26, 27, 21, 26, 27, 28, 22, 27, 28, 29, 23, 28, 29, 14, 30, 31, 33, 24, 30, 31, 32, 34, 25, 31, 32, 35, 30, 33, 34, 36, 31, 33, 34, 35, 37, 32, 34, 35, 38, 33, 36, 37, 43, 34, 36, 37, 38, 44, 35, 37, 38, 39, 45, 38, 39, 40, 45, 39, 40, 41, 45, 40, 41, 42, 45, 41, 42, 46, 36, 43, 44, 37, 43, 44, 45, 38, 39, 40, 41, 44, 45, 46, 42, 45, 46];
Gnzval = [1.3333333333333333, -0.6666666666666666, -0.6666666666666666, -0.6666666666666666, 1.0666666666666667, -0.4, -0.4, 1.0857142857142859, -0.2857142857142857, -0.4, -0.2857142857142857, 0.8412698412698412, -0.2222222222222222, -0.3333333333333333, -0.2222222222222222, 0.5079365079365079, -0.2857142857142857, -0.6666666666666666, 1.0666666666666667, -0.4, -0.4, 1.5666666666666667, -0.5, -0.6666666666666666, -0.3333333333333333, -0.5, 2.0, -0.5, -0.6666666666666666, -0.2857142857142857, -0.5, 1.4523809523809523, -0.6666666666666666, -0.4, 0.6857142857142857, -0.2857142857142857, -0.6666666666666666, 2.6666666666666665, -1.0, -1.0, -0.6666666666666666, -1.0, 3.6666666666666665, -1.0, -1.0, -0.6666666666666666, -1.0, 2.6666666666666665, -1.0, -0.2857142857142857, 1.4228937728937727, -0.4, -0.3333333333333333, -0.25, -0.15384615384615385, -0.4, 2.4, -1.0, -1.0, -1.0, -1.0, 3.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, 3.0, -1.0, -0.3333333333333333, -1.0, 2.0, -0.6666666666666666, -1.0, 2.6666666666666665, -1.0, -0.6666666666666666, -1.0, -1.0, 2.9166666666666665, -0.25, -0.6666666666666666, -0.25, 0.7222222222222222, -0.25, -0.2222222222222222, -0.25, 0.9166666666666666, -0.6666666666666666, -0.25, -0.6666666666666666, 2.083333333333333, -0.5, -0.6666666666666666, -0.5, 1.6666666666666665, -0.5, -0.6666666666666666, -0.6666666666666666, -0.5, 1.6666666666666665, -0.5, -0.6666666666666666, -0.5, 1.6666666666666665, -0.5, -0.2222222222222222, -0.5, 1.2222222222222223, -0.5, -0.6666666666666666, -0.5, 1.1666666666666665, -0.15384615384615385, 0.5371794871794872, -0.25, -0.13333333333333333, -0.6666666666666666, -0.25, 2.9166666666666665, -1.0, -1.0, -0.6666666666666666, -1.0, 2.6666666666666665, -1.0, -0.13333333333333333, 0.47320261437908495, -0.2222222222222222, -0.11764705882352941, -1.0, -0.2222222222222222, 3.2222222222222223, -1.0, -1.0, -1.0, -1.0, 3.0, -1.0, -0.11764705882352941, 0.42291021671826623, -0.2, -0.10526315789473684, -1.0, -0.2, 3.2, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, 2.6666666666666665, -1.0, -0.6666666666666666, -1.0, 2.5, -1.0, -0.5, -1.0, 2.4, -1.0, -0.4, -1.0, 1.3333333333333333, -0.3333333333333333, -0.10526315789473684, 0.28708133971291866, -0.18181818181818182, -1.0, -0.18181818181818182, 2.1818181818181817, -1.0, -1.0, -0.6666666666666666, -0.5, -0.4, -1.0, 3.7888888888888888, -0.2222222222222222, -0.3333333333333333, -0.2222222222222222, 0.5555555555555556];
G = SparseMatrixCSC(Gm, Gn, Gcolptr, Growval, Gnzval);


text:
"""
The **pcl_octree** library provides efficient methods for creating a
hierarchical tree data structure from point cloud data.

http://docs.pointclouds.org/trunk/group__octree.html

## Exports

$(EXPORTS)
"""
module PCLOctree

using DocStringExtensions
using LibPCL
using PCLCommon
using Cxx

const libpcl_octree = LibPCL.find_library_e("libpcl_octree")
try
    Libdl.dlopen(libpcl_octree, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

end # module

ext: jl


text:
function squareseigenvalues(n,m,k)
    if m < 0
        throw(DomainError(m,"order must be non-negative"));
    end
    if n < 1
        throw(DomainError(n,"dimension must be positive"));
    end
    if k > m
        throw(DomainError(k,"degree must be lesser than or equal to m"));
    end
    # for i=0, the eigenvalue is always 0
    values = [1.0];
    # i>1
    for i in 1:1:k
        # might require using log to calculate in case the values get too extreme
        lambda = (2*(m-i)/(2*m+2*i+n))*values[i];
        if m==i
            lambda = factorial(m)/prod([(j+(2*m+n)/2) for j in 1:1:m]);
        end
        push!(values,lambda);
    end
    #normalization = 2*pi^(n/2)/SpecialFunctions.gamma(n/2)
    return values #.* normalization
end

function squaresevaluation(m,p,vars)
    desc = harmonicdecomposition(p,vars)
    n = length(vars)
    k = floor(Int,maxdegree(p)/2)
    if 2*k != maxdegree(p)
        throw(DomainError(p,"polynomial p be of even degree"))
    end
    if 2*k != mindegree(p)
        throw(DomainError(p,"polynomial p must be homogeneous"))
    end
    eigvals = squareseigenvalues(n,m,k)
    norm = sum(vars .* vars)
    return sum( eigvals[floor(Int,u[1]/2)+1]*norm^(floor(Int,u[1]/2))*u[2] for u in desc )
end

function squaresinverseeval(m,p,vars)
    desc = harmonicdecomposition(p,vars)
    n = length(vars)
    k = floor(Int,maxdegree(p)/2)
    if 2*k != maxdegree(p)
        throw(DomainError(p,"polynomial p be of even degree"))
    end
    if 2*k != mindegree(p)
        throw(DomainError(p,"polynomial p must be homogeneous"))
    end
    eigvals = squareseigenvalues(n,m,k)
    norm = sum(vars .* vars)
    return sum( (1/eigvals[floor(Int,u[1]/2)+1])*norm^(floor(Int,u[1]/2))*u[2] for u in desc )
end


text:
# wraps variables found in THCALC.INC for use in Julia

struct THCALC_s
    CUSTOMINTERVALS::Array{Int32,1}
    INTERVALS::OffsetArrays.OffsetArray{Int32,1,Array{Int32,1}}
    function THCALC_s()
        as_lthet = cglobal((:as_lthet_, libaswing), Int32);
        CUSTOMINTERVALS  = unsafe_wrap(Array, as_lthet, 1)

        as_ithet = cglobal((:as_ithet_, libaswing), Int32);
        intervals  = unsafe_wrap(Array, as_ithet, (NBRKX+2))
        INTERVALS = OffsetArrays.OffsetArray(intervals, 0:NBRKX+1)

        return new(CUSTOMINTERVALS, INTERVALS)
    end
end


text:
#Examples:

"""
    cbioproc_example1()

A 3D example without land mask or MeshArrays.
"""
function cbioproc_example1()
    cbioproc_example1("")
end

function cbioproc_example1(dirIn::String)
    in=fill(1.,(90,1170,50))
    SPM,lon,lat=read_SPM(dirIn)
    siz=size(lon)
    out=MatrixInterp(in,SPM,siz)
    #
    outLoop,outName=loop_exampleB(1:13,SPM,siz)
    tmp=transpose(outLoop[1,1])
    heatmap(vec(lon[:,1]),vec(lat[1,:]),tmp)
    return outLoop,outName
end

"""
    cbioproc_example2()

A 2D example with MeshArrays + mask or loop.
"""
function cbioproc_example2()
    dirIn=""
    cbioproc_example2(dirIn)
end

function cbioproc_example2(dirIn::String)
    #
    GCMGridSpec()
    GCMGridLoad()
    msk2d=mask(view(MeshArrays.hFacC,:,:,1),NaN,0)
    SPM,lon,lat=read_SPM(dirIn)
    siz=size(lon)
    #
    in=read_bin(dirIn*"diags/state_2d_set1.0000000732.data",Float32)
    in=convert2gcmfaces(in[:,:,1].*msk2d);
    outMsk=MatrixInterp(in,SPM,siz)
    heatmap(vec(lon[:,1]),vec(lat[1,:]),transpose(outMsk))
    #
    outLoop=loop_exampleA(1:13,SPM,siz);
    heatmap(vec(lon[:,1]),vec(lat[1,:]),transpose(outLoop[1]))
    #
    return outLoop,outMsk
end

"""
    loop_exampleA(indx,SPM,siz)

Applies `MatrixInterp` in a loop over a subset of model output files (`filList[indx]`)
"""
function loop_exampleA(indx,SPM,siz)
   dirIn="diags/"
   filIn="state_2d_set1"
   tmp1=readdir(dirIn)
   tmp1=filter(x -> occursin(filIn,x),tmp1)
   filList=filter(x -> occursin(".data",x),tmp1)
   filList=filList[indx]
   nf=length(filList)
   #show(filList)

   out=Array{Any,1}(undef,nf)
   for ff=1:nf
      fil=dirIn*filList[ff]
      fid = open(fil);
      #println(fil)
      fld = Array{Float32,2}(undef,(90,1170));
      #fld = Array{prec,1}(undef,(n1*n2*n3));
      read!(fid,fld);
      fld = hton.(fld);
      #println(size(fld))
      out[ff]=MatrixInterp(fld,SPM,siz)
   end

   return out
end

"""
    loop_exampleB(indx,SPM,siz)

Applies `MatrixInterp` in a loop over a subset of model output files (`filList[indx]`)
"""
function loop_exampleB(indx,SPM,siz)
   dirIn="diags/"
   filIn="state_2d_set1"
   tmp1=readdir(dirIn)
   tmp1=filter(x -> occursin(filIn,x),tmp1)
   filList=filter(x -> occursin(".data",x),tmp1)
   filList=filList[indx]
   nf=length(filList)
   #show(filList)

   MetaFile=dirIn*filList[1]
   MetaFile=MetaFile[1:end-5]*".meta"
   MetaFile=MetaFileRead(MetaFile)
   nv=length(MetaFile["fldList"])
   nd=MetaFile["nDims"]
   dims=Int.(MetaFile["dimList"][:,1])
   prec=MetaFile["dataprec"]

   out=Array{Any,2}(undef,(nf,nv))
   tmp = Array{prec,nd}(undef,Tuple(dims));
   for ff=1:nf
      fil=dirIn*filList[ff]
      fid = open(fil);
      #println(fil)
      #tmp = Array{prec,1}(undef,(n1*n2*n3));
      for vv=1:nv
         read!(fid,tmp);
         tmp = hton.(tmp);
         out[ff,vv]=MatrixInterp(tmp,SPM,siz)
      end
   end

   return out,MetaFile["fldList"]
end

"""
    prep_MTRX()

Repackage interpolation matrix, mask, etc to `.jld` file.
"""
function prep_MTRX()
    GCMGridSpec()
    GCMGridLoad()
    msk2d=mask(view(MeshArrays.hFacC,:,:,1),NaN,0)
    msk3d=mask(MeshArrays.hFacC,NaN,0)
    msk2d=convert2gcmfaces(msk2d)
    msk3d=convert2gcmfaces(msk3d)

    dirIn=""
    MTRX,lon,lat=read_SPM(dirIn)
    lon=vec(lon[:,1])
    lat=vec(lat[1,:])

    fid = open("GRID_LLC90/RC.data")
    dep=Array{Float64,1}(undef,50)
    read!(fid,dep)
    close(fid)
    dep = -hton.(dep)

    siz2d=(length(lon),length(lat))
    siz3d=(length(lon),length(lat),50)

    save(dirIn*"MTRX.jld", "MTRX", MTRX, "lon", lon, "lat", lat, "dep", dep,
    "msk2d", msk2d, "msk3d", msk3d, "siz2d", siz2d, "siz3d", siz3d)
end

"""
    read_SPM(dirIn::String)

Reads pre-computed interpolation (sparse matrix) from
`dirIn*"interp_precomputed.mat"`.
"""
function read_SPM(dirIn::String)
    #vars = matread(dirIn*"interp_precomputed.mat")
    file = matopen(dirIn*"interp_precomputed.mat")
    interp=read(file, "interp")
    lon=read(file, "lon")
    lat=read(file, "lat")
    SPM=interp["SPM"]
    #println(keys(interp))
    close(file)
    return SPM,lon,lat
end

max_issues_repo_name: darwinproject/CBIOMES-Processing-legacy

max_forks_repo_name: darwinproject/CBIOMES-Processing-legacy


text:
function initialize(M,N,η)
    Ntot=M+N
    Nx=N
    state=EPState{Float64}(Ntot,Nx)
    if η<1 # noisy scenario
       	@extract state A y Σ v av va a μ b s
        b[1:Nx] .*= 1/Nx
        state=EPState{Float64}(A,y,Σ,v,av,va,a,μ,b,s)
    end
    return state
end

function normalize_weights(w,ret,N)
    teacher = w./norm(w);
    student = ret.av[1:N]./norm(ret.av[1:N]);
    return teacher, student
end

function signflip(o,η)
    M=length(o); println("M=$M");
    K_flip=floor(Int64,(1-η)*M); #number of flipped labels
    println("K_flip=$K_flip");
    #flipped = 1:K_flip; o[flipped]*=-1
    flipped = randperm(M)[1:K_flip]; o[flipped]*=-1
    #for i=1:M
    #    if rand() > η
    #        o[i] *= -1
    #    end
    #end
    return o
end

function Lambda_prior(η)
    @assert 0.0<=η<=1.0
    if η==1.0
        println("Noiseless case")
        prior="theta"
    else
        println("Noisy case")
        prior="theta-mixture"
    end
    return prior
end

function Reporter(w)
    function report(iter,state,Δav,Δva,epsconv,maxiter,H,P0)
        av=state.av
        teacher=w./norm(w)
        N=length(w)
        student=av[1:N]./norm(av[1:N])
        overlap=teacher'student
        MSE=Distances.msd(student,teacher)
        MSEdB=10.0*log10(MSE)
        if iter % 100 == 0
            println("$iter $Δav $MSEdB $overlap")
        end
        Δav < epsconv && return true
    end
end

function EP_Regress_ZeroTemp(ξ,w,η;λ=1.0,ρ0=0.25,λ0=1.0,δρ=0.0,δλ=0.0,
η0=0.5,δη=0.0,damp=0.9,maxiter=10000,epsconv=1.0e-4,offset=0.0,callback=(x...)->nothing,
maxvar=1e50,minvar=1e-50,EPInputState::Union{EPState{Float64},Nothing}=nothing,
prior::String="theta_mixture")
    K = length(findall(w.!=0.0))
    M,N = size(ξ)
    o = sign.(ξ*w)
    o_new = signflip(o,η)
    ξ.*=o_new
    Γ = SpikeSlabPrior(ρ0,λ0,δρ,δλ)
    (prior=="theta") && (Λ=ThetaPrior())
    (prior=="theta_mixture") && (Λ=ThetaMixturePrior(η0,δη))
    P0 = [ [Γ for i=1:N]; [Λ for i=1:M] ]
    if typeof(EPInputState) == Nothing
        EPInputState = EPState{Float64}(M+N,N)
    end
    callback==Reporter(w) && println("Printing (iter, Δav, MSEdB, overlap):")
    runtime =
        @elapsed ret=
        expectation_propagation(Vector{Term{Float64}}(undef,0),
        P0,ξ,offset*ones(M),damp=damp,maxiter=maxiter,epsconv=epsconv,maxvar=maxvar,
        minvar=minvar,callback=callback,state=EPInputState,inverter=x->inv(cholesky(Hermitian(x))))
    teacher,student=normalize_weights(w,ret,N);
    Overlap=teacher'student; println("α=$(M/N), $(ret.converged)");
    MSEdB=10.0*log10(Distances.msd(teacher,student))
    return ret, (Γ, Λ), runtime, MSEdB, Overlap
end

max_stars_repo_name: mpieropan/EPDilutedClassifier

max_issues_repo_name: mpieropan/EPDilutedClassifier

max_forks_repo_name: mpieropan/EPDilutedClassifier


text:
#julia %USERPROFILE%\.julia\packages\PackageCompiler\oT98U\juliac.jl -C ivybridge -vRe gbk2UTF8Bom.jl
using StringEncodings

function TransferGbkCSV2UTF8BomCsv( ;srcFile::String ,TgrFile::String )
    f = open(srcFile, enc"GB2312", "r");
    fout = open(TgrFile, enc"UTF-8", "w");
    write(fout,UInt8(239) )
    write(fout,UInt8(187) )
    write(fout,UInt8(191) )
    for line in readlines( f )
        println( line )
        write( fout,line )
        write(fout,"\r\n")
    end
    close( fout)
    close(f)
end
Base.@ccallable function julia_main1()::Cint
    srcFile = "D:/Jl-2019-0521/Jl/glbCfg/xxxxx-ParasGBk-FULL.csv"
    TgrFile  = "GYYS-001-AlarmParas.csv"
    TransferGbkCSV2UTF8BomCsv( srcFile= srcFile ,TgrFile= TgrFile )
    return 0
end
#julia_main1()

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PredicatePartition(pred)

A method for partitioning spatial objects with a given predicate
function `pred(i, j)`.
"""
struct PredicatePartition <: PredicatePartitionMethod
  pred::Function
end

(p::PredicatePartition)(i, j) = p.pred(i, j)


text:
snoop_file = joinpath(@__DIR__, "test-snoop.txt")
ENV["PkgCompilerEnabled"] = "1"

using PkgCompiler
PkgCompiler.set_global_snoop_file(snoop_file)
PkgCompiler.set_pkg_compiler_dir(@__DIR__)

include("code.jl")

PkgCompiler.build_system_image()


text:
export
    Proposal,
    bootstrapprop

struct Proposal
    mu0::Union{Float,Vector{Float}}
    noise::Function
    mean::Function
    loglik::Function
end

function bootstrapprop(g::GaussianHMM, mu0::Union{Float,Vector{Float}}=0.0)
    hmm = HMM(g)
    n = nothing
    Proposal(
        (mu0==0.0) ? (g.dimx>1?zeros(g.dimx): mu0 ) : mu0,
        (g.dimx>1)?(k=n)->(g.cholQ'*randn(g.dimx)):((k=n)->g.cholQ'*randn()),
        (k=n,xkm1=n,yk=n)      -> hmm.transmean(k, xkm1),
        (k=n,xkm1=n,yk=n,xk=n) -> hmm.transloglik(k, xkm1, xk)
    )
end

# # reverese bootstrap (reversed dynamic)
# function rbootstrap(lg::LinearGaussian)
#     hmm = HMM(lg)
#     n   = nothing
#     pia = pinv(lg.A)
#     Proposal(
#         (k=n)                  -> pia * lg.cholQ' * randn(lg.dimx)
#         (k=n,xkp1=n,yk=n)      -> pia * xkp1
#         (k=n,xk=n,yk=n,xkp1)   -> hmm.transloglik(k, xk, xkp1)
#     )
# end


text:
function incremental(f,min::Int64,max::Int64,N::Int64)
    t = range(min,stop = max,length = N) #makes the domain
    intervalo = Array{Float64}(undef,N-1,2) #prepares all intervals
    intervalo[:,1], intervalo[:,2] = t[1:end-1], t[2:end]
    return intervalo[(f.(t[1:end-1])) .* (f.(t[2:end])) .< 0 ,:] #chooses all intervals that the product of the edges is negative
end

function incremental(f,min::Float64,max::Float64,N::Int64)
    t = range(min,stop = max,length = N) #makes the domain
    intervalo = Array{Float64}(undef,N-1,2) #prepares all intervals
    intervalo[:,1], intervalo[:,2] = t[1:end-1], t[2:end]
    return intervalo[(f.(t[1:end-1])) .* (f.(t[2:end])) .< 0 ,:]
end


text:
module MortalityTables
using OffsetArrays
using Parsers
using QuadGK
using Requires
using UnPack
using XMLDict
using Pkg.Artifacts

include("table_source_map.jl")
include("MetaData.jl")
include("death_distribution.jl")
include("MortalityTable.jl")
include("XTbML.jl")
include("get_SOA_table.jl")
include("parameterized_models.jl")
include("life_expectancy.jl")

table_dirs = Dict(
    "mort.soa.org" => artifact"mort.soa.org",
)

export MortalityTable,
    survival,
    decrement,
    life_expectancy,
    omega,
    TableMetaData,
    SelectMortality,
    UltimateMortality,
    Balducci,
    Uniform,
    Constant,
    DeathDistribution,
    get_SOA_table,
    Makeham, Gompertz, MakehamGompertz,
    hazard,cumhazard,
    mortality_vector

# lazy load part of the package
function __init__()
    @require CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" include("CSV.jl")
end

end # module


text:
## ---------  Transfer regression

function transfer_regr(ds :: DataSettings)
    # Regression deviations use RegressionTables as inputs
    fPath = data_file(raw_transfer_regr(ds));
    rt = read_regression_file(fPath);
    rename_regressors(rt);
    return rt
end

## -----------  Tuition regression (net price)
# on quality, parental, gpa

function tuition_regr(ds :: DataSettings)
    # Regression deviations use RegressionTables as inputs
    fPath = data_file(raw_tuition_regr(ds));
    rt = read_regression_file(fPath);
    rename_regressors(rt);
    # Adjust all means and std errors to match model units
    # modelUnits  &&  convert_to_model_dollars!(rt);
    return rt
end

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


text:
### I Am here. move this to NMRSpectraSimulator.

using LinearAlgebra, FFTW
import BSON, Statistics, Random
import PyPlot

import NMRSpectraSimulator

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

# for loading something with Interpolations.jl
import OffsetArrays
import Interpolations

import PlotlyJS
import Plots
Plots.plotly()

import Destruct

include("./helpers/final_helpers.jl")
include("./helpers/resonance_helpers.jl")

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

#projects_dir = save_folder_path
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/L-Serine-700"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/D-(+)-Glucose-NRC-600"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/D-(+)-Glucose-700"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/L-Isoleucine-700"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/L-Leucine-500-2mM"
projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/L-Glutamine-700"

plots_save_folder = joinpath(projects_dir, "plots")
isdir(plots_save_folder) || mkdir(plots_save_folder)
project_title = "resonance groupings"

T = Float64

### load block.
#load_path = joinpath(joinpath(projects_dir, project_name), "results_full.bson")
load_path = joinpath(projects_dir, "results_full.bson")
dict = BSON.load(load_path)
As = collect( dict[:As][i] for i = 1:length(dict[:As]) )
Es = collect( NMRSpectraSimulator.κCompoundFIDType(As[i]) for i = 1:length(As) )

Δsys_cs = convert(Vector{Vector{Float64}}, dict[:Δsys_cs])
y = convert(Vector{Complex{Float64}}, dict[:y])
U_y = convert(Vector{Float64}, dict[:U_y])
SW = dict[:SW]
fs = dict[:fs]
ν_0ppm = dict[:ν_0ppm]

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)
P_y = hz2ppmfunc.(U_y)

ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(As, hz2ppmfunc)
ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))
u_offset = 0.5
u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)

p_star_set = dict[:p_star_set]
κ_lb_default = dict[:κ_lb_default]
κ_ub_default = dict[:κ_ub_default]
κ_star_set = dict[:κ_star_set]
d_star_set = dict[:d_star_set]
β_star_set = dict[:β_star_set]
λ_star_set = dict[:λ_star_set]
cost_inds_set = dict[:cost_inds_set]
w = dict[:w]

As = collect( Es[n].core for n = 1:length(Es))
### end load block.

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

### plot.

f = uu->NMRSpectraSimulator.evalmixture(uu, As)

# test params.
ΩS_ppm = collect( hz2ppmfunc.( NMRSpectraSimulator.combinevectors(A.Ωs) ./ (2*π) ) for A in As )
#ΩS_ppm_flat = NMRSpectraSimulator.combinevectors(ΩS_ppm)

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)

## parameters that affect qs.
# A.ss_params.d, A.ss_params.κs_λ, A.ss_params.κs_β
# A.d_singlets, A.αs_singlets, A.Ωs_singlets, A.β_singlets, A.λ0, A.κs_λ_singlets
q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As[1:1])

#g = uu->evalitpproxycompoundresonance(uu, As[1])
A = As[1]
g1 = uu->evalitpproxysysresonance(A.qs, uu, A.ss_params.d, A.ss_params.κs_λ, A.ss_params.κs_β)
g2 = uu->evalsingletsresonance(uu, A.d_singlets, A.αs_singlets, A.Ωs_singlets,
A.β_singlets, A.λ0, A.κs_λ_singlets)

f_U = f.(U)
q_U = q.(U)
g1_U = g1.(U)
g2_U = g2.(U)

#g_qs_U, g_singlets_U = convertresonancetimeseries(g_U)

z1 = Destruct.destruct(g1_U);
q_singlets_evals = Destruct.destruct(g2_U);

# collect( length(z1[i]) for i = 1:length(z1))
# collect( length(z2[i]) for i = 1:length(z2))
qs_evals = assembleqsevals(z1, A.qs)

sanity_check = zeros(Complex{Float64}, length(q_U))
if !isempty(z1)
    sanity_check += sum(z1)
end
if !isempty(q_singlets_evals)
    sanity_check += sum(q_singlets_evals)
end
println("sanity: ", norm(sanity_check - q_U))
@assert norm(sanity_check - q_U) < 1e-14
println()

discrepancy = abs.(f_U-q_U)
max_val, ind = findmax(discrepancy)
println("relative discrepancy = ", norm(discrepancy)/norm(f_U))
println("max discrepancy: ", max_val)
println()

## visualize.
PyPlot.plot(P, real.(f_U), label = "f")
PyPlot.plot(P, real.(q_U), label = "q")

#visualizesubgroups(P, qs_evals, q_singlets_evals, real)

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("f vs q")

#@assert 1==2

plots_save_path = joinpath(plots_save_folder, "groups_real.html")
title_string = "$(project_title), real"
savefiggroups(plots_save_path, title_string,
    P, q_U, qs_evals, q_singlets_evals, real)

#
plots_save_path = joinpath(plots_save_folder, "groups_imag.html")
title_string = "$(project_title), imag"
savefiggroups(plots_save_path, title_string,
    P, q_U, qs_evals, q_singlets_evals, real)
#
plots_save_path = joinpath(plots_save_folder, "groups_modulus.html")
title_string = "$(project_title), modulus"
savefiggroups(plots_save_path, title_string,
    P, q_U, qs_evals, q_singlets_evals, real)


text:
export anp_1d

"""
    anp_1d(;
        dim_embedding::Integer,
        num_encoder_layers::Integer,
        num_encoder_heads::Integer,
        num_decoder_layers::Integer,
        noise_type::String="het",
        pooling_type::String="mean",
        σ::Float32=1f-2,
        learn_σ::Bool=true
    )

# Arguments
- `dim_embedding::Integer`: Dimensionality of the embedding.
- `num_encoder_layers::Integer`: Number of layers in the encoder.
- `num_encoder_heads::Integer`: Number of heads in the encoder.
- `num_decoder_layers::Integer`: Number of layers in the decoder.
- `noise_type::String="het"`: Type of noise model. Must be "fixed", "amortised", or "het".
- `pooling_type::String="mean"`: Type of pooling. Must be "mean" or "sum".
- `σ::Float32=1f-2`: Initialisation of the fixed observation noise.
- `learn_σ::Bool=true`: Learn the fixed observation noise.

# Returns
- `Model`: Corresponding model.
"""
function anp_1d(;
    dim_embedding::Integer,
    num_encoder_layers::Integer,
    num_encoder_heads::Integer,
    num_decoder_layers::Integer,
    noise_type::String="het",
    pooling_type::String="mean",
    σ::Float32=1f-2,
    learn_σ::Bool=true
)
    dim_x = 1
    dim_y = 1
    num_noise_channels, noise = build_noise_model(
        dim_y       =dim_y,
        noise_type  =noise_type,
        pooling_type=pooling_type,
        σ           =σ,
        learn_σ     =learn_σ
    )
    return Model(
        Parallel(
            Chain(
                InputsCoder(),
                DeterministicLikelihood()
            ),
            Chain(
                attention(
                    dim_x             =dim_x,
                    dim_y             =dim_y,
                    dim_embedding     =dim_embedding,
                    num_heads         =num_encoder_heads,
                    num_encoder_layers=num_encoder_layers
                ),
                DeterministicLikelihood()
            ),
            Chain(
                MLPCoder(
                    batched_mlp(
                        dim_in    =dim_x + dim_y,
                        dim_hidden=dim_embedding,
                        dim_out   =dim_embedding,
                        num_layers=num_encoder_layers
                    ),
                    batched_mlp(
                        dim_in    =dim_embedding,
                        dim_hidden=dim_embedding,
                        dim_out   =2dim_embedding,
                        num_layers=num_encoder_layers
                    )
                ),
                HeterogeneousGaussianLikelihood()
            )
        ),
        Chain(
            Materialise(),
            batched_mlp(
                dim_in    =dim_x + 2dim_embedding,
                dim_hidden=dim_embedding,
                dim_out   =num_noise_channels,
                num_layers=num_decoder_layers
            ),
            noise
        )
    )
end


text:
export RoyaltyModelNoHet, RoyaltyModel, ObservationRoyalty, DataRoyalty,
    observed_choices

# Model structure
#---------------------------

"Royalty rates"
struct RoyaltyModelNoHet <: AbstractRoyaltyModel end
struct RoyaltyModel      <: AbstractRoyaltyModel end

# Abstract data strucutres
#---------------------------

struct ObservationRoyalty{M<:AbstractRoyaltyModel,I<:Integer, T<:Real, V<:AbstractVector{T}} <: AbstractObservation
    model::M
    y::I
    x::V
    xbeta::T
    num_choices::I
    function ObservationRoyalty(model::M,y::I,x::V,xbeta::T,num_choices::I) where {M<:AbstractRoyaltyModel,I<:Integer, T<:Real, V<:AbstractVector{T}}
        isfinite(xbeta) || throw(DomainError(xbeta, "$xbeta = xβ = $x * β not finite")) # FIXME
        0 < y <= num_choices || throw(DomainError(y, "y = $y should be 1 <= y <= $num_choices"))
        return new{M,I,T,V}(model, y, x, xbeta, num_choices)
    end
end

struct DataRoyalty{M<:AbstractRoyaltyModel,I<:Integer,AV<:AbstractVector, T<:Real} <: AbstractDataSet
    model::M
    y::Vector{I}
    x::Matrix{T}
    xbeta::Vector{T}
    choices::AV
    function DataRoyalty(model::M, y::Vector{I}, x::Matrix{T},choices::AV) where {M,I,AV,T}
        k,n = size(x)
        length(y) == n  || throw(DimensionMismatch())
        issorted(choices) || throw(error("choices not sorted"))
        L = length(choices)
        extrema(y) == (1,L) || throw(error("lowest royalty ID > 1"))
        L == length(countmap(y)) || throw(error("highest royalty ID < L = $L"))
        return new{M,I,AV,T}(model, y, x, Vector{T}(undef, n), choices)
    end
end

const DataOrObsRoyalty = Union{ObservationRoyalty{M},DataRoyalty{M}} where {M}
const ObservationGroupRoyalty = ObservationGroup{<:DataRoyalty}
const AbstractDataStructureRoyalty = Union{ObservationRoyalty,DataRoyalty,ObservationGroupRoyalty}

function DataRoyalty(m::AbstractRoyaltyModel, d::DataRoyalty)
    DataRoyalty(m, _y(d), _x(d), choices(d))
end

# Common interfaces
#---------------------------

_xbeta(      d::DataOrObsRoyalty) = d.xbeta
choices(d::DataRoyalty) = d.choices
num_choices(d::ObservationRoyalty) = d.num_choices
num_choices(d::DataRoyalty) = length(choices(d))

observed_choices(d::DataRoyalty) = getindex(choices(d), _y(d))

@deprecate _num_choices(d::DataOrObsRoyalty) num_choices(d)

# DataRoyalty interface
#---------------------------

length(   d::DataRoyalty) = length(_y(d))
group_ptr(d::DataRoyalty) = OneTo(length(d)+1)
obs_ptr(  d::DataRoyalty) = OneTo(length(d)+1)

function choice_in_model(d::DataOrObsRoyalty, l::Integer)
    0 < l <= num_choices(d) && return true
    throw(DomainError(l, "$l outside of 1:$(num_choices(d))"))
end

function choice_in_model(d, ls)
    for l in ls
        choice_in_model(d,l)
    end
    return nothing
end

# Iteration over data
#---------------------------

function Observation(d::DataRoyalty, k::Integer)
    k ∈ eachindex(d) || throw(BoundsError(d,k))
    y = getindex(_y(d), k)
    x = view(_x(d), :, k)
    xbeta = getindex(_xbeta(d), k)
    return ObservationRoyalty(_model(d), y, x, xbeta, num_choices(d))
end

getindex(   g::ObservationGroupRoyalty, k) = Observation(_data(g), getindex(grouprange(g), k))
Observation(g::ObservationGroupRoyalty, k) = getindex(g,k)

# Observation
#---------------------------

function ==(a::ObservationRoyalty, b::ObservationRoyalty)
    _y(a) == _y(b) &&
    _x(a) == _x(b) &&
    _xbeta(a) == _xbeta(b) &&
    num_choices(a) == num_choices(b)
end

# Updating data
#---------------------------

function update_xbeta!(d::DataRoyalty, theta::AbstractVector)
    # length(theta) == _num_x(d) || throw(DimensionMismatch("theta: $(size(theta)) and x: $(size(_x(d)))"))
    mul!(_xbeta(d), _x(d)', theta)
    return nothing
end

update!(d::DataRoyalty, theta) = update_xbeta!(d,theta_royalty_β(d,theta))

# Parameter interfaces
#---------------------------

# length of RoyaltyModel parameters & gradient
extra_parm(::RoyaltyModel) = 2
extra_parm(::RoyaltyModelNoHet) = 0
extra_parm(d::DataOrObsRoyalty{<:AbstractRoyaltyModel}) = extra_parm(_model(d))

_nparm(d::DataOrObsRoyalty) = _num_x(d) + num_choices(d) - 1 + extra_parm(d)

idx_royalty(d::Union{DataOrObsRoyalty,AbstractRoyaltyModel}) = OneTo(_nparm(d))
theta_royalty(d, theta) = view(theta, idx_royalty(d))

# parameter vector
idx_royalty_ρ(d::Union{DataOrObsRoyalty{<:RoyaltyModelNoHet},RoyaltyModelNoHet}) = 1:0
idx_royalty_ψ(d::Union{DataOrObsRoyalty{<:RoyaltyModelNoHet},RoyaltyModelNoHet}) = 1:0
idx_royalty_β(d::DataOrObsRoyalty{<:RoyaltyModelNoHet}) = (1:_num_x(d))
idx_royalty_κ(d::DataOrObsRoyalty{<:RoyaltyModelNoHet}) = _num_x(d) .+ (1:num_choices(d)-1)

function idx_royalty_κ(d::DataOrObsRoyalty{RoyaltyModelNoHet}, l)
    choice_in_model(d,l)
    return _num_x(d) .+ l
end

idx_royalty_ρ(d::Union{DataOrObsRoyalty{RoyaltyModel},RoyaltyModel}) = 1
idx_royalty_ψ(d::Union{DataOrObsRoyalty{RoyaltyModel},RoyaltyModel}) = 2
idx_royalty_β(d::DataOrObsRoyalty{RoyaltyModel}) = 2 .+ (1:_num_x(d))
idx_royalty_κ(d::DataOrObsRoyalty{RoyaltyModel}) = 2 + _num_x(d) .+ (1:num_choices(d)-1)
function idx_royalty_κ(d::DataOrObsRoyalty{RoyaltyModel}, l::Integer)
    choice_in_model(d,l)
    return 2 + _num_x(d) + l
end
function idx_royalty_κ(d::DataOrObsRoyalty{RoyaltyModel}, ls::AbstractVector)
    choice_in_model(d,ls)
    return (2 + _num_x(d)) .+ ls
end

# get coefs
theta_royalty_ρ(d, theta) = theta[idx_royalty_ρ(d)]
theta_royalty_ψ(d, theta) = theta[idx_royalty_ψ(d)]
theta_royalty_β(d, theta) = view(theta, idx_royalty_β(d))
theta_royalty_κ(d, theta) = view(theta, idx_royalty_κ(d))
theta_royalty_κ(d, theta, l::Integer) = theta[idx_royalty_κ(d,l)]
theta_royalty_κ(d, theta, ls) = view(theta, idx_royalty_κ(d,ls))

kappa_level_to_cumsum(x) = vcat(first(x), sqrt.(2 .* diff(x)) )
kappa_cumsum_to_level(x) = first(x) .+ vcat(0, 0.5 .* cumsum(x[2:end].^2))
function update_kappa_level_to_cumsum!(x)
    x .= kappa_level_to_cumsum(x)
end

function theta_royalty_level_to_cumsum(d,x)
    xnew = copy(x)
    kap = theta_royalty_κ(d,x)
    xnew[idx_royalty_κ(d)] .= kappa_level_to_cumsum(kap)
    return xnew
end
function theta_royalty_cumsum_to_level(d,x)
    xnew = copy(x)
    kap = theta_royalty_κ(d,x)
    xnew[idx_royalty_κ(d)] .= kappa_cumsum_to_level(kap)
    return xnew
end

# @deprecate theta_roy(d, theta) theta_royalty(d,theta)

# check if theta is okay
theta_royalty_check(d, theta) = issorted(theta_royalty_κ(d,theta))

function check_theta(d::DataOrObsRoyalty, theta)
    _nparm(d) == length(theta) || throw(DimensionMismatch("_nparm(d) = $(_nparm(d)) != length(theta) = $(length(theta))"))
    all(isfinite.(theta)) || throw(error("theta = $theta not finite"))
    kappa = theta_royalty_κ(d,theta)
    # issorted(kappa) || @warn "kappa = $kappa not sorted"
    return true
end

function coefnames(d::DataRoyalty)
    nms = Vector{String}(undef, _nparm(d))
    if _model(d) isa RoyaltyModel
        nms[idx_royalty_ρ(d)] = "\\rho"
        nms[idx_royalty_ψ(d)] = "\\psi^0"
    end
    for (i,nmi) in enumerate(idx_royalty_β(d))
        nms[nmi] = "\\beta_$i"
    end
    for (i, nmi) in enumerate(idx_royalty_κ(d))
        nms[nmi] = "\\kappa_$i"
    end
    return nms
end

"""
    DataRoyalty(u,v,X,theta,L)

Simulate dataset for `RoyaltyModel` using `u,v` to make `ψ1`
"""
function DataRoyalty(u::AbstractVector, v::AbstractVector, X::Matrix, theta::Vector, rates::AbstractVector)

    k,nobs = size(X)
    L = length(rates)
    L >= 3 || throw(error("L = $L !>= 3"))

    k == length(theta) - (L-1) - 2 || throw(DimensionMismatch("dim mismatch b/w X: $(size(X)) and θ = $theta"))
    k >= 1 || throw(error("theta too short"))
    nobs == length(u) == length(v) || throw(DimensionMismatch())

    # get ψ1
    ψ1 = similar(u)
    dψ1dρ = similar(u)
    update_ψ1!(ψ1, u, v, first(theta))
    update_dψ1dθρ!(dψ1dρ, u, v, first(theta))

    eps    = randn(nobs)

    rstar  = theta[2] .* ψ1 .+ X'*theta[2 .+ (1:k)] .+ eps
    cutoffs = kappa_cumsum_to_level(theta[end-L+2:end])
    l = map((r) ->  searchsortedfirst(cutoffs, r), rstar)
    data = DataRoyalty(RoyaltyModel(),l,X,rates)

    return data
end

function DataRoyalty(u, v, X, theta, L::Integer=3)
    if L == 3
        rates = [1/8, 3/16, 1/4]
    else
        throw(error("L != 3"))
    end
    return DataRoyalty(u,v,X,theta,rates)
end

"""
    DataRoyalty(u,v,theta,L)

Simulate dataset for `RoyaltyModel` generating random X
"""
function DataRoyalty(u::AbstractVector, v::AbstractVector, theta::Vector, L::Integer=3)
    k = length(theta) - (L-1) - 2
    nobs = length(u)
    X = randn(k, nobs)
    return DataRoyalty(u,v,X,theta,L)
end

function DataRoyalty(num_i::Integer, theta::Vector, L::Integer=3)
    return DataRoyalty(randn(num_i), randn(num_i), theta, L)
end


text:
# This file was generated, do not modify it. # hide
RidgeRegressor = @load RidgeRegressor pkg="MultivariateStats"
KNNRegressor = @load KNNRegressor

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
###################
# OpSum/DualOpSum #
###################
abstract type AbsOpSum{P,N,T} <: DiracOp{P,N} end

const OpDict{N,T} = Dict{OpLabel{N},T}

mutable struct OpSum{P,N,T} <: AbsOpSum{P,N,T}
    ptype::P
    dict::OpDict{N,T}
    # OpSum(ptype, dict) = new(ptype, dict)
    OpSum(ptype, dict::OpDict{0}) = error("Cannot construct a 0-factor operator; did you mean to construct a scalar?")
end

OpSum(ptype::P, dict::OpDict{N,T}) where {P,N,T} = OpSum{P,N,T}(ptype, dict)
OpSum(kt::Ket{P,N,A}, br::Bra{P,N,B}) where {P,N,A,B} = OpSum(ptype(kt), cons_outer!(OpDict{N,promote_type(A,B)}(), kt, br))

function cons_outer!(result, kt, br)
    for (k,kc) in dict(kt)
        for (b,bc) in dict(br)
            newc = kc * bc'
            if newc != 0
                result[OpLabel(k, b)] = newc
            end
        end
    end
    return result
end

struct DualOpSum{P,N,T} <: AbsOpSum{P,N,T}
    op::OpSum{P,N,T}
end

# DualOpSum(op::OpSum{P,N,T}) where {P,N,T} = DualOpSum{P,N,T}(op)
DualOpSum(items...) = DualOpSum(OpSum(items...))

Base.convert(::Type{OpSum}, opc::DualOpSum) = eager_ctran(opc.op)
Base.promote_rule(::Type{O}, ::Type{D}) where {O<:OpSum, D<:DualOpSum} = OpSum

######################
# Accessor functions #
######################
dict(op::OpSum) = op.dict
dict(opc::DualOpSum) = dict(opc.op)

ptype(op::OpSum) = op.ptype
ptype(opc::DualOpSum) = ptype(opc.op)

#######################
# Dict-Like Functions #
#######################
Base.eltype(::AbsOpSum{P,N,T}) where {P,N,T} = T

Base.copy(op::OpSum) = OpSum(ptype(op), copy(dict(op)))
Base.copy(opc::DualOpSum) = DualOpSum(copy(opc.op))

Base.similar(op::OpSum, d=similar(dict(op)); P=ptype(op)) = OpSum(P, d)
Base.similar(opc::DualOpSum, d=similar(dict(opc)); P=ptype(opc)) = DualOpSum(P, d)

Base.:(==)(a::OpSum{P,N}, b::OpSum{P,N}) where {P,N} = ptype(a) == ptype(b) && dict(filternz(a)) == dict(filternz(b))
Base.:(==)(a::DualOpSum{P,N}, b::DualOpSum{P,N}) where {P,N} = a.op == b.op
Base.:(==)(a::DiracOp, b::DiracOp) = ==(promote(a,b)...)

Base.hash(op::AbsOpSum) = hash(dict(filternz(op)), hash(ptype(op)))
Base.hash(op::AbsOpSum, h::UInt64) = hash(hash(op), h)

Base.length(op::AbsOpSum) = length(dict(op))

Base.getindex(op::OpSum, label::OpLabel) = op.dict[label]
Base.getindex(op::OpSum, k::StateLabel, b::StateLabel) = op.dict[OpLabel(k,b)]
Base.getindex(opc::DualOpSum, label::OpLabel) = opc.op[label']'
Base.getindex(opc::DualOpSum, k::StateLabel, b::StateLabel) = opc.op[OpLabel(b,k)]'
Base.getindex(op::AbsOpSum, k, b) = op[StateLabel(k), StateLabel(b)]

Base.setindex!(op::OpSum, c, label::OpLabel) = (op.dict[label] = c)
Base.setindex!(op::OpSum, c, k::StateLabel, b::StateLabel) = (op.dict[OpLabel(k,b)] = c)
Base.setindex!(opc::DualOpSum, c, label::OpLabel) = (opc.op[label'] = c')
Base.setindex!(opc::DualOpSum, c, k::StateLabel, b::StateLabel) = (opc.op[OpLabel(b,k)] = c')
Base.setindex!(op::AbsOpSum, c, k, b) = setindex!(op, c, StateLabel(k), StateLabel(b))

Base.haskey(op::OpSum, label::OpLabel) = haskey(dict(op), label)
Base.haskey(opc::DualOpSum, label::OpLabel) = haskey(opc.op, label')
Base.haskey(op::AbsOpSum, k, b) = haskey(op, OpLabel(k, b))

Base.get(op::OpSum, label::OpLabel, default=predict_zero(eltype(op))) = get(dict(op), label, default)
Base.get(opc::DualOpSum, label::OpLabel, default=predict_zero(eltype(opc))) = get(dict(opc), label, default')'
Base.get(op::AbsOpSum, k, b, default=predict_zero(eltype(op))) = get(op, OpLabel(k, b), default)

Base.delete!(op::OpSum, label::OpLabel) = (delete!(dict(op), label); return op)
Base.delete!(opc::DualOpSum, label::OpLabel) = delete!(opc.op, label')
Base.delete!(op::AbsOpSum, k, b) = delete!(op, OpLabel(k, b))

Base.collect(op::OpSum) = collect(dict(op))
Base.collect(opc::DualOpSum{P,N,T}) where {P,N,T} = collect_pairs!(Array(@compat(Tuple{OpLabel{N}, T}), length(opc)), opc)

function collect_pairs!(result, opc::DualOpSum)
    i = 1
    for (k,v) in dict(opc)
        result[i] = ctpair(k,v)
        i += 1
    end
    return result
end

Base.iterate(op::AbsOpSum) = iterate(dict(op))
Base.iterate(op::OpSum, i) = iterate(dict(op), i)

function Base.iterate(opc::DualOpSum, i)
    (k,v), n = iterate(dict(opc), i)
    return ((k',v'), n)
end

# Base.done(op::AbsOpSum, i) = done(dict(op), i)
Base.first(op::AbsOpSum) = iterate(op, iterate(op))

#############
# ctranpose #
#############
eager_ctran(op::OpSum) = similar(op, Dict(collect(op')))

Base.adjoint(op::OpSum) = DualOpSum(op)
Base.adjoint(opc::DualOpSum) = opc.op

#########
# inner #
#########
function inner(br::Bra{P,N,A}, op::OpSum{P,N,B}) where {P,N,A,B}
    prodtype = ptype(op)
    result = StateDict{N, inner_coefftype(br, op)}()
    return Bra(prodtype, inner_load!(result, br, op, prodtype))
end

function inner(op::OpSum{P,N,A}, kt::Ket{P,N,B}) where {P,N,A,B}
    prodtype = ptype(op)
    result = StateDict{N, inner_coefftype(op, kt)}()
    return Ket(prodtype, inner_load!(result, op, kt, prodtype))
end

function inner(a::OpSum{P,N,A}, b::OpSum{P,N,B}) where {P,N,A,B}
    prodtype = ptype(a)
    result = OpDict{N, inner_coefftype(a, b)}()
    return OpSum(prodtype, inner_load!(result, a, b, prodtype))
end

function inner(a::OpSum{P,N,A}, b::DualOpSum{P,N,B}) where {P,N,A,B}
    prodtype = ptype(a)
    result = OpDict{N, inner_coefftype(a, b)}()
    return OpSum(prodtype, inner_load!(result, a, b, prodtype))
end

function inner(a::DualOpSum{P,N,A}, b::OpSum{P,N,B}) where {P,N,A,B}
    prodtype = ptype(a)
    result = OpDict{N, inner_coefftype(a, b)}()
    return OpSum(prodtype, inner_load!(result, a, b, prodtype))
end

inner(br::Bra, opc::DualOpSum) = inner(opc.op, br')'
inner(opc::DualOpSum, kt::Ket) = inner(kt', opc.op)'
inner(a::DualOpSum, b::DualOpSum) = inner(b.op, a.op)'

function inner_load!(result, a::OpSum, b::OpSum, prodtype)
    for (o1,v) in dict(a), (o2,c) in dict(b)
        add_to_dict!(result,
                     OpLabel(klabel(o1), blabel(o2)),
                     inner_mul(v, c, prodtype, blabel(o1), klabel(o2)))
    end
    return result
end

function inner_load!(result, a::OpSum, b::DualOpSum, prodtype)
    for (o1,v) in dict(a), (o2,c) in dict(b)
        add_to_dict!(result,
                     OpLabel(klabel(o1), klabel(o2)),
                     inner_mul(v, c', prodtype, blabel(o1), blabel(o2)))
    end
    return result
end

function inner_load!(result, a::DualOpSum, b::OpSum, prodtype)
    for (o1,v) in dict(a), (o2,c) in dict(b)
        add_to_dict!(result,
                     OpLabel(blabel(o1), blabel(o2)),
                     inner_mul(v', c, prodtype, klabel(o1), klabel(o2)))
    end
    return result
end

function inner_load!(result::Dict{K,T}, br::Bra, op::OpSum, prodtype) where {K,T}
    for (o,v) in dict(op)
        add_to_dict!(result, blabel(o), brcoeff(dict(br), prodtype, klabel(o), v, T))
    end
    return result
end

function inner_load!(result::Dict{K,T}, op::OpSum, kt::Ket, prodtype) where {K,T}
    for (o,v) in dict(op)
        add_to_dict!(result, klabel(o), ktcoeff(dict(kt), prodtype, blabel(o), v, T))
    end
    return result
end

function brcoeff(brdict, prodtype, klabel, v, ::Type{T}) where {T}
    coeff = predict_zero(T)
    for (blabel,c) in brdict
        coeff += inner_mul(c', v, prodtype, klabel, blabel)
    end
    return coeff'
end

function ktcoeff(ktdict, prodtype, blabel, v, ::Type{T}) where {T}
    coeff = predict_zero(T)
    for (klabel,c) in ktdict
        coeff += inner_mul(c, v, prodtype, klabel, blabel)
    end
    return coeff
end

Base.:*(br::Bra, op::DiracOp) = inner(br,op)
Base.:*(op::DiracOp, kt::Ket) = inner(op,kt)
Base.:*(a::DiracOp, b::DiracOp) = inner(a,b)

###################################
# Functional Operator Application #
###################################
struct DualFunc
    f::Function
end

Base.:*(op::Function, kt::Ket) = op(kt)
Base.:*(br::Bra, op::Function) = op(br)
Base.:*(op::DualFunc, kt::Ket) = (kt' * op.f)'
Base.:*(br::Bra, op::DualFunc) = (op.f * br')'

Base.adjoint(f::Function) = DualFunc(f)
Base.adjoint(fc::DualFunc) = fc.f

##############
# act/act_on #
##############
act_on(op::AbsOpSum, br::Bra, i) = act_on(op', br', i)'

# clear up ambiguity warnings
act_on(op::OpSum{P,1}, kt::Ket{P,1}, i) where P = i==1 ? inner(op, kt) : throw(BoundsError())
act_on(opc::DualOpSum{P,1}, kt::Ket{P,1}, i) where P = i==1 ? inner(opc, kt) : throw(BoundsError())

function act_on(op::OpSum{P,1,A}, kt::Ket{P,N,B}, i) where {P,N,A,B}
    prodtype = ptype(op)
    result = StateDict{N, inner_coefftype(op, kt)}()
    return Ket(prodtype, act_on_dict!(result, op, kt, i, prodtype))
end

function act_on(op::DualOpSum{P,1,A}, kt::Ket{P,N,B}, i) where {P,N,A,B}
    prodtype = ptype(op)
    result = StateDict{N, inner_coefftype(op, kt)}()
    return Ket(prodtype, act_on_dict!(result, op, kt, i, prodtype))
end

function act_on_dict!(result, op::OpSum, kt::Ket, i, prodtype)
    for (o,c) in dict(op), (k,v) in dict(kt)
        add_to_dict!(result,
                     setindex(k, klabel(o)[1], i),
                     inner_mul(c, v, prodtype, blabel(o)[1], k[i]))
    end
    return result
end

function act_on_dict!(result, op::DualOpSum, kt::Ket, i, prodtype)
    for (o,c) in dict(op), (k,v) in dict(kt)
        add_to_dict!(result,
                     setindex(k, blabel(o)[1], i),
                     inner_mul(c', v, prodtype, klabel(o)[1], k[i]))
    end
    return result
end

##########
# tensor #
##########
tensor(a::OpSum{P}, b::OpSum{P}) where {P} = OpSum(ptype(a), tensordict(dict(a), dict(b)))
tensor(a::DualOpSum, b::DualOpSum) = tensor(a.opc, b.opc)'
tensor(a::DiracOp, b::DiracOp) = tensor(promote(a,b)...)

Base.:*(kt::Ket, br::Bra) = tensor(kt,br)

###########
# Scaling #
###########
scale!(op::OpSum, c::Number) = (dscale!(dict(op), c); return op)
scale!(c::Number, op::OpSum) = scale!(op, c)
scale!(opc::DualOpSum, c::Number) = DualOpSum(scale!(opc.op, c'))
scale!(c::Number, opc::DualOpSum) = scale!(opc, c)

scale(op::OpSum, c::Number) = similar(op, dscale(dict(op), c))
scale(c::Number, op::OpSum) = scale(op, c)
scale(opc::DualOpSum, c::Number) = DualOpSum(scale(opc.op, c'))
scale(c::Number, opc::DualOpSum) = scale(opc, c)

Base.:*(c::Number, op::DiracOp) = scale(c, op)
Base.:*(op::DiracOp, c::Number) = scale(op, c)
Base.:/(op::DiracOp, c::Number) = scale(op, 1/c)

###########
# + and - #
###########
Base.:-(op::OpSum) = scale(-1, op)
Base.:-(opc::DualOpSum) = DualOpSum(-opc.op)

Base.:+(a::OpSum{P,N}, b::OpSum{P,N}) where {P,N} = similar(b, add_merge(dict(a), dict(b)))
Base.:-(a::OpSum{P,N}, b::OpSum{P,N}) where {P,N} = similar(b, sub_merge(dict(a), dict(b)))

Base.:+(a::DualOpSum{P,N}, b::DualOpSum{P,N}) where {P,N} = DualOpSum(a.op + b.op)
Base.:-(a::DualOpSum{P,N}, b::DualOpSum{P,N}) where {P,N} = DualOpSum(a.op - b.op)

Base.:+(a::DiracOp, b::DiracOp) = +(promote(a,b)...)
Base.:-(a::DiracOp, b::DiracOp) = a + (-b)

#################
# Normalization #
#################
norm(op::OpSum) = sqrt(sum(abs2, values(dict(op))))
norm(opc::DualOpSum) = norm(opc.op)

normalize(op::DiracOp) = scale(1/norm(op), op)
normalize!(op::DiracOp) = scale!(1/norm(op), op)

#########
# Trace #
#########
function LinearAlgebra.tr(op::OpSum)
    result = predict_zero(eltype(op))
    for (o,v) in dict(op)
        if klabel(o)==blabel(o)
            result += v
        end
    end
    return result
end

LinearAlgebra.tr(opc::DualOpSum) = tr(opc.op)'

#################
# Partial trace #
#################
ptrace(op::DiracOp{P,1}, over) where {P} = over == 1 ? tr(op) : throw(BoundsError())
ptrace(op::DiracOp{P,N}, over) where {P,N} = OpSum(ptype(op), ptrace_dict!(OpDict{N-1,eltype(op)}(), op, over))

function ptrace_dict!(result, op::OpSum, over)
    for (o,v) in dict(op)
        if klabel(o)[over] == blabel(o)[over]
            add_to_dict!(result, traceout(o, over), v)
        end
    end
    return result
end

function ptrace_dict!(result, opc::DualOpSum, over)
    for (o,v) in dict(opc)
        if blabel(o)[over] == klabel(o)[over]
            add_to_dict!(result, traceout_dual(o, over), v')
        end
    end
    return result
end

#####################
# Partial Transpose #
#####################
ptranspose(op::DiracOp{P,N}, over) where {P,N} = OpSum(ptype(op), ptrans_dict!(OpDict{N,eltype(op)}(), op, over))

function ptrans_dict!(result, op::OpSum, over)
    for (o,v) in dict(op)
        add_to_dict!(result, ptranspose(o, over), v)
    end
    return result
end

function ptrans_dict!(result, opc::DualOpSum, over)
    for (o,v) in dict(opc)
        add_to_dict!(result, ptranspose_dual(o, over), v')
    end
    return result
end

########################
# Misc. Math Functions #
########################
nfactors(::AbsOpSum{P,N}) where {P,N} = N
xsubspace(op::AbsOpSum, x) = similar(op, filter((k,v)->is_sum_x(k,x), dict(op)))
filternz!(op::AbsOpSum) = (filter!(nzcoeff, dict(op)); return op)
filternz(op::AbsOpSum) = similar(op, filter(nzcoeff, dict(op)))
switch(op::AbsOpSum, i, j) = maplabels(label->switch(label,i,j), op)
permute(op::AbsOpSum, perm::Vector) = maplabels(label->permute(label,perm), op)

purity(op::DiracOp) = tr(op^2)
commute(a::DiracOp, b::DiracOp) = (a*b) - (b*a)
anticommute(a::DiracOp, b::DiracOp) = (a*b) + (b*a)

inner_eval(f, op::DiracOp) = mapcoeffs(x->inner_eval(f,x),op)

function matrep(op::DiracOp, labels)
    T = promote_type(inner_rettype(ptype(op)), eltype(op))
    return T[bra(i) * op * ket(j) for i in labels, j in labels]
end

function matrep(op::DiracOp, labels...)
    iter = product(labels...)
    T = promote_type(inner_rettype(ptype(op)), eltype(op))
    return T[bra(i...) * op * ket(j...) for i in iter, j in iter]
end

function matrep(op::Union{DualFunc, Function}, labels)
    return [bra(i) * op * ket(j) for i in labels, j in labels]
end

function matrep(op::Union{DualFunc, Function}, labels...)
    iter = Iterators.product(labels...)
    return [bra(i...) * op * ket(j...) for i in iter, j in iter]
end

######################
# Printing Functions #
######################
labelrepr(op::OpSum, o::OpLabel, pad) = "$pad$(op[o]) $(ktstr(klabel(o)))$(brstr(blabel(o)))"
labelrepr(opc::DualOpSum, o::OpLabel, pad) = "$pad$(opc[o']) $(ktstr(blabel(o)))$(brstr(klabel(o)))"

Base.summary(op::DiracOp) = "$(typeof(op)) with $(length(op)) operator(s)"
Base.show(io::IO, op::AbsOpSum) = dirac_show(io, op)
# Base.showcompact(io::IO, op::AbsOpSum) = dirac_showcompact(io, op)
Base.repr(op::AbsOpSum) = dirac_repr(op)

export OpSum,
    ptrace,
    ptranspose,
    matrep,
    xsubspace,
    nfactors,
    maplabels!,
    mapcoeffs!,
    mapcoeffs,
    maplabels,
    filternz,
    filternz!,
    purity,
    act_on,
    inner_eval


text:
export RQ, precode!, ltgenerate

"""

Return the parameter tuple `(Kp, J, S, H, W)` for RaptorQ codes with `K` source symbols.
"""
function RQ_parameters(K::Integer)
    0 < K <= 56403 || throw(DomainError(K, "K must be in [1, 56403]"))
    i = searchsortedfirst(view(RQ_parameter_table, :, 1), K)
    Kp = RQ_parameter_table[i, 1]
    J = RQ_parameter_table[i, 2]
    S = RQ_parameter_table[i, 3]
    H = RQ_parameter_table[i, 4]
    W = RQ_parameter_table[i, 5]
    Kp, J, S, H, W
end

"""
    RQ

RaptorQ code.

"""
struct RQ <: AbstractErasureCode
    K::Int # number of source symbols
    Kp::Int # number of source symbols including padding
    J::Int # systematic index
    S::Int # number of LDPC symbols
    H::Int # number of HDPC symbols
    W::Int # number of LT symbols
    L::Int # =Kp+S+H
    P::Int # number of PI symbols
    P1::Int # smallest prime >= P
    U::Int # number of PI symbols that are not HDPC symbols
    B::Int # number of LT symbols excluding the LDPC symbols
    function RQ(K::Integer)
        Kp, J, S, H, W = RQ_parameters(K)
        L = Kp + S + H
        P = L - W
        P1 = nextprime(P)
        U = P - H
        B = W - S
        new(K, Kp, J, S, H, W, L, P, P1, U, B)
    end
end
Base.show(io::IO, code::RQ) = print(io, "RQ($(code.K))")

"""

RaptorQ standardized rand function. Maps the integers `y` and `i` to a psuedo-randomly generated
number between `0` and `m-1` (inclusive).
"""
function RQ_rand(y::Integer, i::Integer, m::Integer)
    0 <= y || throw(DomainError(y, "y must be non-negative"))
    0 <= i < 256 || throw(DomainError(i, "i must be in [0, 256)"))
    0 <= m || throw(DomainError(m, "m must be positive"))
    x0 = (y + i) % 256
    x1 = (Int(floor(y / 256)) + i) % 256
    x2 = (Int(floor(y / 65536)) + i) % 256
    x3 = (Int(floor(y / 16777216)) + i) % 256
    result = xor(V0[x0+1], V1[x1+1])
    result = xor(result, V2[x2+1])
    result = xor(result, V3[x3+1])
    result % m
end

"""

RaptorQ degree distribution. Maps an integer `v` in `[0, 2^20)` to a degree, where `v` is assumed
to be uniformly distributed over the range.
"""
function RQ_deg(v::Integer)
    0 <= v < 1048576 || throw(DomainError(v, "v must be in [0, 2^20)"))
    searchsortedlast(RQ_DEGREE_TABLE, v)
end

"""

RaptorQ standardized tuple function. Takes as input the number of source symbols `Kp` and an ESI
`X`, and returns the tuple `(d, a, b, d1, a1, b1)`, which uniquely determines the LT symbol with
ESI `X`.
"""
function RQ_tuple(X::Integer, c::RQ)
    A = 53591 + 997c.J
    if (A % 2 == 0) A += 1 end
    B = 10267*(c.J+1)
    y = (B + X*A) % 2^32
    v = RQ_rand(y, 0, 2^20)
    d = min(RQ_deg(v), c.W-2)
    a = 1 + RQ_rand(y, 1, c.W-1)
    b = RQ_rand(y, 2, c.W)
    if (d < 4) d1 = 2 + RQ_rand(X, 3, 2)
    else d1 = 2 end
    a1 = 1 + RQ_rand(X, 4, c.P1-1)
    b1 = RQ_rand(X, 5, c.P1)
    d, a, b, d1, a1, b1
end

function ldpc_constraint_matrix(c::RQ)
    indices = [Vector{Int}() for _ in 1:c.S]
    for i in 0:c.B-1
        a = 1 + Int(floor(i/c.S))
        b = i % c.S
        push!(indices[b+1], i+1)
        b = (b + a) % c.S
        push!(indices[b+1], i+1)
        b = (b + a) % c.S
        push!(indices[b+1], i+1)
    end
    for i in 0:c.S-1
        a = i % c.P
        b = (i + 1) % c.P
        if c.W+a == i+c.Kp
            push!(indices[i+1], c.W+b+1)
        elseif c.W+b == i+c.Kp
            push!(indices[i+1], c.W+a+1)
        else
            push!(indices[i+1], i+c.Kp+1)
            push!(indices[i+1], c.W+a+1)
            push!(indices[i+1], c.W+b+1)
        end
    end
    hcat([sparsevec(Is, ones(GF256, length(Is)), c.L) for Is in indices]...)
end

function hdpc_constraint(c::RQ, A::Matrix{GF256}, i::Integer)
    Is = Vector{Int}()
    Vs = Vector{GF256}()
    for j in 1:c.Kp+c.S
        if !iszero(A[i, j])
            push!(Is, j)
            push!(Vs, A[i, j])
        end
    end
    push!(Is, c.Kp+c.S+i)
    push!(Vs, one(GF256))
    sparsevec(Is, Vs, c.L)
end

function hdpc_constraint_matrix(c::RQ)
    alpha = GF256(0x02)
    MT = zeros(GF256, c.H, c.Kp+c.S)
    GAMMA = zeros(GF256, c.Kp+c.S, c.Kp+c.S)
    for j in 0:c.Kp+c.S-2
        i = RQ_rand(j+1, 6, c.H)
        MT[i+1,j+1] = one(GF256)
        i = (i + RQ_rand(j+1, 7, c.H-1) + 1) % c.H
        MT[i+1,j+1] = one(GF256)
    end
    j = c.Kp+c.S
    for i in 1:c.H
        MT[i,j] = alpha^(i-1)
    end
    for i in 0:c.Kp+c.S-1
        for j in 0:i
            GAMMA[i+1,j+1] = alpha^(i-j)
        end
    end
    A::Matrix{GF256} = MT*GAMMA
    return hcat([hdpc_constraint(c, A, i) for i in 1:c.H]...)

    # TODO: more efficient implementation. exploits the structure of the matrix.
    # (currently doesn't work)
    # A = zeros(GF256, c.H, c.Kp+c.S)
    # for j in 0:c.Kp+c.S-2 # col index of MT
    #     i = RQ_rand(j+1, 6, c.H) # row index of MT
    #     l = j # row index of GAMMA
    #     for k in 0:l # col index of GAMMA
    #         A[i+1,k+1] += alpha^(l-k)
    #         N[i0+i][2][k+1] += alpha^(l-k)
    #     end
    #     i = (i + RQ_rand(j+1, 7, c.H-1) + 1) % c.H
    #     for k in 0:l
    #         A[i+1,k+1] += alpha^(l-k)
    #         N[i0+i][2][k+1] += alpha^(l-k)
    #     end
    # end
    # j = c.S+c.H-1
    # l = j
    # for i in 0:c.H-1
    #     for k in 0:i-1
    #         A[i+1,k+1] += alpha^(i+l-k)
    #         N[i0+i][2][k+1] += alpha^(i+l-k)
    #     end
    # end
end

"""

Return the constraint matrix composed of the LDPC and HDPC constraint.
"""
precode_constraint_matrix(code::RQ) = hcat(ldpc_constraint_matrix(code), hdpc_constraint_matrix(code))

"""

Return a `SparseVector` corresponding to the LT symbol with ESI `X`.
"""
function lt_constraint(c::RQ, X::Integer)
    d, a, b, d1, a1, b1 = RQ_tuple(X, c)
    Is = zeros(Int, d+d1)
    Is[1] = b+1
    for j in 1:d-1
        b = (b + a) % c.W
        Is[j+1] = b+1
    end
    while b1 >= c.P
        b1 = (b1+a1) % c.P1
    end
    Is[d+1] = c.W+b1+1
    for j in 1:d1-1
        b1 = (b1 + a1) % c.P1
        while b1 >= c.P
            b1 = (b1+a1) % c.P1
        end
        Is[d+j+1] = c.W+b1+1
    end
    sparsevec(Is, ones(GF256, d+d1), c.L)
end

generator_matrix(code::RQ, Xs::AbstractVector{<:Integer}) = hcat([lt_constraint(code, X) for X in Xs]...)
constraint_matrix(code::RQ, Xs::AbstractVector{<:Integer}) = hcat(precode_constraint_matrix(code), generator_matrix(code, Xs))

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
using MySQL
using DataFrames

conn = Nothing

if length(ARGS) !=
    println("Usage:\n\tjulia -L cli.jl <host> <username> <password> <port> <database>")
    println("""\nOR\n    conn = MySQL.connect("127.0.0.1", "dalla", "password", port = 3306, db = "dalla_stats")""")
else
    conn = MySQL.connect(ARGS[1], ARGS[2], ARGS[3], port = ARGS[4], db = ARGS[5])
end

include("DallaStats.jl")

println("Usage:\n\tDallaStats.update_device_table(conn, month::Integer)\n\n")
# DallaStats.update_device_table(conn, 11)

max_stars_repo_name: egeldenhuys/dalla-stats-cli

max_issues_repo_name: egeldenhuys/dalla-stats-cli

max_forks_repo_name: egeldenhuys/dalla-stats-cli


text:
module ReverseDiffTape

importall Base

# package code goes here
function __init__()
    # Logging.configure(level=DEBUG)
    @show "loading ReverseDiffTape.jl"
end

export @timing, @asserting,
#constant
    OP, S_TO_OC,
#types
    TYPE_V, TYPE_P, TYPE_O,
#Objects
    AD, AD_O, AD_P, AD_V, Tape, HessStorage,
#Functions
    append_array,
    tapeBuilder, tapeBuilderNoHess,appendMultParam,appendTapeMultParam,buildSumTape, #building tape from Julia expression
    mergeTapes,getMaxWorkingSize,resizeHessStorage,
    feval,
    grad_reverse, grad_reverse_dense, grad_structure,
    hess_structure, hess_reverse, reset_hess,
    tape_report_mem

@inline function append_array{I,V}(dest::Vector{V},d_offset::I,src::Vector{V},s_offset::I, n::I)
    for i=1:n
        @inbounds dest[i+d_offset] = src[i+s_offset]
    end
end

include("./types.jl")
include("./operator.jl")
include("./func_eval.jl")
include("./reverse_grad.jl")
include("./reverse_hess_ep.jl")

end # module


text:
#klubovi.jl
using Printf

struct Klub
    ime::String
    brojPoena::Int
    datiGolovi::Int
    primljeniGolovi::Int
end

klubovi = []
push!(klubovi, Klub("Benevento", 30, 30, 55))
push!(klubovi, Klub("Milan", 63, 57, 35))
push!(klubovi, Klub("Juventus", 62, 61, 27))
push!(klubovi, Klub("Napoli", 59, 65, 34))
push!(klubovi, Klub("Inter", 74, 69, 27))
push!(klubovi, Klub("Roma", 54, 54, 44))
push!(klubovi, Klub("Fiorentina", 30, 38, 49))
push!(klubovi, Klub("Lazio", 55, 46, 38))
push!(klubovi, Klub("Atalanta", 61, 71, 38))

function uporedi(klub1, klub2)
	return klub1.brojPoena > klub2.brojPoena || klub1.brojPoena == klub2.brojPoena &&
	       (klub1.datiGolovi-klub1.primljeniGolovi>klub2.datiGolovi-klub2.primljeniGolovi)
end

function ispisKlubova(klubovi)
  for i in 1:length(klubovi)
     @printf("%2i. %-15s %3i %3i %3i\n", i, klubovi[i].ime,
              klubovi[i].brojPoena, klubovi[i].datiGolovi, klubovi[i].primljeniGolovi)
  end
end


text:
# Julia wrapper for header: /usr/include/libavutil/rc4.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_rc4_alloc,
    av_rc4_init,
    av_rc4_crypt

function av_rc4_alloc()
    ccall((:av_rc4_alloc, libavutil), Ptr{AVRC4}, ())
end

function av_rc4_init(d, key, key_bits::Integer, decrypt::Integer)
    ccall((:av_rc4_init, libavutil), Cint, (Ptr{AVRC4}, Ptr{UInt8}, Cint, Cint), d, key, key_bits, decrypt)
end

function av_rc4_crypt(d, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_rc4_crypt, libavutil), Cvoid, (Ptr{AVRC4}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), d, dst, src, count, iv, decrypt)
end


text:
# load diabetes.data

using LinearAlgebra
using CSV
using DataFrames
using Plots

Diabete_Data = CSV.read("diabetes.csv",DataFrame)

n,d = size(Diabete_Data)
d = d-1

#println(Diabete_Data)

# compute the optimal theta

X = zeros(n,d)
y = zeros(n,1)

for i=1:d
    X[:,i] = Diabete_Data[:,i]
end
#X = [ones(n,1)  X[:,3]]
X = [ones(n,1) X]
y = Diabete_Data[:,end]

theta_opt = X\y

# achieved loss
MSE = norm(X*theta_opt-y,2)^2/n

println("theta_opt: ", theta_opt)
println("RMSE: ",sqrt(MSE))

# plot
figure()
plot(y,y,"k")
plot(y,X*theta_opt,"o",alpha=0.4)
xlabel(L"$y$")
ylabel(L"predicted $y$")
title("prediction accuracy")
axis("square")
grid("on")

figure()
hist(y-X*theta_opt, bins=10, density=true)
title("prediction error")
grid("on")

# put JHK data into your predictor

#age
#sex
#bmi
#map	   mean arterial pressure
#s1  tc :   total cholesterol
#s2  ldl    low density lipoprotein
#s3  hdl    high density lipoprotein
#s4  tch
#s5  ltg
#s6  glu
#        age  sex  bmi    map   tc     ldl    hdl     tch   ltg   glu
X_JHK = [41   1    18.3   90    171    80.0   74.9    2     4.75  90.0]
X_JHK = [1 X_JHK]

y_JHK = X_JHK*theta_opt

# plot

figure()
plot(y,y,"k")
plot(y,X*theta_opt,"o",alpha=0.4)
plot(y_JHK,y_JHK,"ro",alpha=0.4)
xlabel(L"$y$")
ylabel(L"predicted $y$")
axis("square")
grid("on")

max_stars_repo_path: Primer Predavanje/diabetes.jl

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_path: Primer Predavanje/diabetes.jl

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_path: Primer Predavanje/diabetes.jl

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using LinearAlgebra, Interpolations, Plots
using HA_solver

#SOme simple tests, not even a proper tests, just checking if function works

#initialize model
Model = HA_solver.Ayiagari_params()

#define paramters:

#grid definition: number of points
Model.a_n = 500
Model.ϵ_n = 9
Model.a_min = 0.0 #minimal number of points
#expotential grid: params for construction
Model.Δ = 0.01
Model.g = 0.015
#definition of the income process
Model.rho = 0.6
Model.sigma_e = 0.4
Model.MC = "T" #use Tauchen method
#preferences paramters
Model.β = 0.96
Model.σ = 5.0
#production function
Model.δ = 0.08
Model.α = 0.36
#solution updating and initail value of the interest rate
Model.γ = 0.95
Model.start_r =0.02

σ_set = [1.0, 3.0, 5.0]
ρ_set = [0.0, 0.3, 0.6, 0.9]
sigma_earnings = [0.4]
    #orginal solution from Ayiagari
    println(eachindex(σ_set))
r_computed = zeros(3,4,2)
for s in eachindex(sigma_earnings)
    for i in eachindex(σ_set)

        for j in eachindex(ρ_set)
            Model.σ = σ_set[i]
            Model.rho = ρ_set[j]
            if(i ==1)
                Model.γ = 0.99
                Model.start_r =  0.04
                if(j ==3)
                    Model.start_r =  0.035
                end
            else
                Model.γ = 0.97
                Model.start_r = 0.02
                if(i==3 & s==2)
                    Model.start_r = 0.0
                    Model.γ = 0.95
                end
            end

            Model.sigma_e = sigma_earnings[s]
            println(i)
            println(σ_set[i])
            r_computed[i,j,s] = HA_solver.Find_eq_Ayiagari(Model)

        end
    end
end

write("replicated_sigma=02", r_computed[:,:,1])
write("replicated_sigma=04", r_computed[:,:,2])


text:
include("guiding_center_4d_settings.jl")

set_config(:quasi_newton_refactorize, 1)
set_config(:nls_atol, 2eps())
set_config(:nls_rtol, 2eps())
set_config(:nls_stol, 2eps())
set_config(:nls_nmin, 2)
set_config(:nls_nmax, 20)


text:
@testset "IndicatorDifferentFrom2PseudoMILP: $(fct_type), $(T), $(A)" for fct_type in ["vector of variables", "vector affine function"], T in [Int, Bool, Float64], A in [MOI.ACTIVATE_ON_ZERO, MOI.ACTIVATE_ON_ONE]
    mock = MOIU.MockOptimizer(AbsoluteValueIndicatorPseudoMILPModel{T}())
    model = if A == MOI.ACTIVATE_ON_ZERO
        COIB.Indicator0DifferentFrom2PseudoMILP{T}(mock)
    elseif A == MOI.ACTIVATE_ON_ONE
        COIB.Indicator1DifferentFrom2PseudoMILP{T}(mock)
    end

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        MOI.Indicator{A, MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        MOI.Indicator{A, CP.DifferentFrom{T}},
    )

    x, _ = MOI.add_constrained_variable(model, MOI.ZeroOne())
    if T == Int
        y, _ = MOI.add_constrained_variable(model, MOI.Integer())
    elseif T == Bool
        y, _ = MOI.add_constrained_variable(model, MOI.ZeroOne())
    elseif T == Float64
        y = MOI.add_variable(model)
    else
        @assert false
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x, y])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x, y])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, MOI.Indicator{A}(CP.DifferentFrom(zero(T))))

    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, CP.DifferentFrom{T}}(1)] # Why a 1 for this specific case? Usually, it's -1.

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.DifferentFrom{T}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]

        if T == Int
            @test MOIB.added_constraint_types(typeof(bridge)) == [
                (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
                (MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}}),
            ]
        elseif T == Bool
            @test MOIB.added_constraint_types(typeof(bridge)) == [
                (MOI.VectorAffineFunction{Bool}, MOI.Indicator{A, MOI.EqualTo{Bool}}),
            ]
        elseif T == Float64
            @test MOIB.added_constraint_types(typeof(bridge)) == [
                (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
                (MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}}),
            ]
        else
            @assert false
        end

        @test MOI.get(bridge, MOI.NumberOfVariables()) == ((T == Bool) ? 0 : 1)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.AbsoluteValue}()) == ((T == Bool) ? 0 : 1)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}}}()) == ((T == Float64) ? 1 : 0)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}}}()) == ((T == Int) ? 1 : 0)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.EqualTo{T}}}()) == ((T == Bool) ? 1 : 0)

        if T != Bool
            @test MOI.get(bridge, MOI.ListOfVariableIndices()) == [bridge.var_abs]
        end
        if T != Bool
            @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.AbsoluteValue}()) == [bridge.con_abs]
        end
        if T == Float64
            @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}}}()) == [bridge.con_indic]
        end
        if T == Int
            @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}}}()) == [bridge.con_indic]
        end
        if T == Bool
            @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.EqualTo{T}}}()) == [bridge.con_indic]
        end
    end

    if T != Bool
        @testset "Absolute value" begin
            @test MOI.is_valid(model, bridge.con_abs)
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_abs)
            @test length(f.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.con_abs) == CP.AbsoluteValue()

            t1 = f.terms[1]
            @test t1.output_index == 1
            @test t1.scalar_term.coefficient === one(T)
            @test t1.scalar_term.variable == bridge.var_abs

            t1 = f.terms[2]
            @test t1.output_index == 2
            @test t1.scalar_term.coefficient === one(T)
            @test t1.scalar_term.variable == y
        end
    end

    if T == Float64
        @testset "Strictly greater than" begin
            @test MOI.is_valid(model, bridge.con_indic)
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_indic)
            @test length(f.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.con_indic) == MOI.Indicator{A}(CP.Strictly(MOI.GreaterThan(zero(T))))

            t1 = f.terms[1]
            @test t1.output_index == 1
            @test t1.scalar_term.coefficient == one(T)
            @test t1.scalar_term.variable == x

            t2 = f.terms[2]
            @test t2.output_index == 2
            @test t2.scalar_term.coefficient == one(T)
            @test t2.scalar_term.variable == bridge.var_abs
        end
    end

    if T == Int
        @testset "Greater than" begin
            @test MOI.is_valid(model, bridge.con_indic)
            @test MOI.get(model, MOI.ConstraintSet(), bridge.con_indic) == MOI.Indicator{A}(MOI.GreaterThan(one(T)))
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_indic)
            @test length(f.terms) == 2

            t1 = f.terms[1]
            @test t1.output_index == 1
            @test t1.scalar_term.coefficient == one(T)
            @test t1.scalar_term.variable == x

            t2 = f.terms[2]
            @test t2.output_index == 2
            @test t2.scalar_term.coefficient == one(T)
            @test t2.scalar_term.variable == bridge.var_abs
        end
    end

    if T == Bool
        @testset "EqualTo" begin
            @test MOI.is_valid(model, bridge.con_indic)
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_indic)
            @test length(f.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.con_indic) == MOI.Indicator{A}(MOI.EqualTo(one(T)))

            t1 = f.terms[1]
            @test t1.output_index == 1
            @test t1.scalar_term.coefficient == one(T)
            @test t1.scalar_term.variable == x

            t2 = f.terms[2]
            @test t2.output_index == 2
            @test t2.scalar_term.coefficient == one(T)
            @test t2.scalar_term.variable == y
        end
    end
end


text:
function visualise(t, target=nothing)

   top, right, bottom, left = boundingbox(t, target)

   path = p -> normpath("$(@__DIR__)/$p")
   files = path.(["../assets/js/utils.js", "../assets/js/animate.js", "../assets/css/basic.css"])
   sc = Scope(imports=files)
   c = t.c
   r, ws = c.r, c.w

   bb = Dict(
      "top"=>top,
      "right"=>right,
      "bottom"=>bottom,
      "left"=>left
   )

   fields = Dict(
      "distance"=>[0.0, "m"],
      "height"=>[0.0, "m"],
      "time"=>[0.0, "m"],
      "release_angle"=>[rad2deg(r), "deg"],
      "wind_speed"=>[ws, "m/s"],
   )
   if target != nothing
      fields["target"] = [target, "m"]
   end

   id = sc.id
   onimport(sc,  @js function ()
      createCanvas($(id), "main");
      createOutputBar($(id), "output", $(fields));
      animate($(id), "main", $(t.l), $(t.sol), $(bb), $(target), $(ws))
   end)

   sc(dom"div#_container_"())
end

function boundingbox(t, target)
   col = (y, i) -> map(x -> x[i], y)

   xs = col(t.sol.Projectile, 1)
   ys = col(t.sol.Projectile, 2)

   top = max(ys...)
   right = max(xs...)
   bottom = min(ys...)
   left = min(xs...)
   if target != nothing
      right = max(target, right)
      left = min(target, left)
   end
   return top, right, bottom, left
end


text:
function bt03()
#
# Brugnano and Trigiante
#
    p = poly([im*ones(5);-im*ones(5);0.5im*ones(4);-0.5im*ones(4);0.75im;-0.75im])
    z = [im 5; -im  5; 0.5im  4; -0.5im  4; .75im  1; -0.75im  1];

    p, PolyZeros(z)
end


text:
#!/usr/bin/env julia

using RobotOS
@rosimport barc.msg: ECU,  Z_KinBkMdl
@rosimport data_service.msg: TimeData
@rosimport geometry_msgs.msg: Vector3,  PointStamped
@rosimport geometry_msgs.msg: Vector3
@rosimport ackermann_msgs.msg: AckermannDrive, AckermannDriveStamped
@rosimport nav_msgs.msg: Path, Odometry
@rosimport duckietown_msgs.msg: Twist2DStamped, LanePose

rostypegen()
using barc.msg
using data_service.msg
using geometry_msgs.msg
using ackermann_msgs.msg
using nav_msgs.msg
using JuMP
using Ipopt
using Distances
using duckietown_msgs.msg

# define model parameters
L_a     = 0.1         # distance from CoG to front axel
L_b     = 0.1         # distance from CoG to rear axel
dt      = 0.1           # time step of system
la      =0.08           # center to front wheel of trycicle model
T       = 150  # MPC horizon
dt=0.1 # Smapling time
global Getpos  = false
zmax=[10;10;100*pi]
umax=[0.6;3*dt]
n=3
m=2
# predictive horizon

l       = 0.1   #length of the car
X=[0.0;0.0;0.0]
path = Path()

function solveMPC(l,n,m,T,z0,zT,zmax,umax,dt,step)
    #T : horizon
    #z0 : intial position for every iteration
    #zT : reference   +++++++++++++++
    #n :dim of input
    #m : dim of state
    #
    mpc = Model(solver=IpoptSolver(print_level=0))
    @defVar(mpc, -zmax[i] <= z[t=1:T+1,i=1:n] <= zmax[i])
    @defVar(mpc, -umax[i] <= u[t=1:T+1,i=1:m] <= umax[i])

    # Cost
    #@setObjective(mpc, Min,
    #    sum{100*(z[1,t]+z[2,t]+z[3,t]+z[4,t])^2 + sum{u[j,t]^2,j=1:m},t=0:T})
    #@setObjective(mpc, Min,
    #sum{(z[t,1]-zT[t,1])^2+(z[t,2]-zT[t,2])^2+10e-5*(z[t,3]-zT[t,3])^2  ,t=2:T+1})
@defNLExpr(mpc, bta[i = 1:T], atan( L_a / (L_a + L_b) * tan(u[i,1]) ) )
      @setObjective(mpc, Min,
    sum{(z[T,1]-zT[1])^2+(z[T,2]-zT[2])^2+(z[T,3]-zT[3])^2,t=1:T})
    # Link state and control across the horizon
    for t = 1:T
        @addNLConstraint(mpc, z[t+1,1] == z[t,1] + dt*u[t,2]*cos(z[t,3]))
        @addNLConstraint(mpc, z[t+1,2] == z[t,2] + dt*u[t,2]*sin(z[t,3]))
        @addNLConstraint(mpc, z[t+1,3] == z[t,3] + dt*u[t,2]/l*sin(bta[t]))
        #@addNLConstraint(mpc, z[t+1,4] == z[t,4] + dt*u[t,2])

    end

    # Initial conditions
    @addConstraint(mpc, z[1,:] .== z0)
    # Final state
    #@addConstraint(mpc, z[:,T] .== zT)
    # Solve the NLP
    solve(mpc)
    # Return the control plan
    step = step +1
    return getValue(u), getValue(z), step
end

#------------------- sub state estimation-------------------
function SEcb(msg::Odometry)
    # update mpc initial condition
    global Getpos = true
    X[1] =    msg.pose.pose.position.x
    X[2] =    msg.pose.pose.position.y
    X[3] =    msg.pose.pose.orientation.w
    #X[4] =    msg.twist.twist.linear.x

end

#--------------  Find the nearest point on the reference of the car ----------------
function find_nearest(x,f,start)
    #print(x[1:2])
    #print(f[1:2,1])
    dis = Array(Float64,size(f,1))
    ref = Array(Float64,T)
    #print(size(dis))
    for i in 1:size(f,1)
        dis[i] = euclidean(x[1:2],f[i,1:2])
    end

    refstep=4
    ref_start = indmin(dis)
    if ref_start <= start
       ref_start = start+1
     end
    print(ref_start)

    if (ref_start > size(f,1) - T*refstep) && (ref_start < size(f,1))
        #print(size(f,1)-ref_start)
	r=(size(f,1)-ref_start)%refstep
	if r != 0
                   ref = [ f[ref_start:refstep:size(f,1)-r,:] ; f[refstep-r:refstep:(T*refstep-(size(f,1)-ref_start)+r),:]   ]
               else
	    ref = [ f[ref_start:refstep:size(f,1),:] ; f[1:refstep:1+(T*refstep-(size(f,1)-ref_start)),:]   ]
                end
    elseif ref_start > size(f,1)
        ref_start = ref_start - size(f,1)
        ref = f[ref_start:refstep:ref_start+T*refstep,:]
    elseif ref_start == size(f,1)
	ref = f[1:refstep:T*refstep+1,:]
    else
        ref = f[ref_start:refstep:ref_start+T*refstep,:]
        #print(ref)
    end

    return ref,ref_start
end

#--------------------Main function-----------------------------
function main()
    zT = [0;1.7;1/2*pi]  # Initial state - COO
    Xb =[0;2;pi]
    step=0
    # initiate node, set up publisher / subscriber topic
    init_node("mpc_parking")
    pub_car_cmd = Publisher("mpc/car_cmd", Twist2DStamped, queue_size=1)
    pub_first_refp = Publisher("goal", PointStamped, queue_size=1)
    sub  = Subscriber{Odometry}("euler_odom",  SEcb, queue_size=1)
    c1=1
    loop_rate = Rate(10)

    while ! is_shutdown()
        #print(X)
        if Getpos == true
            # find the nearest point on the path# get reference in horizon
            if step == T
            car_control_msg = Twist2DStamped()
            car_control_msg.v = 0
            car_control_msg.omega = 0
            #print(v)
            #publish commands
            publish(pub_car_cmd, car_control_msg)
                return
            end
            #Xb[1] = X[1] - l/2*cos(pi/2) + l/2*sin(X[3])
            #Xb[2] = X[2] - l/2*sin(pi/2) + l/2*cos(X[3])
            #Xb[3] = X[3]
            ## run mpc, publish command
            u_vec, z_vec, step = solveMPC(l,n,m,T-step,X,zT,zmax,umax,dt, step)

            println("solve!")
            #print(u_vec[1,1] ,u_vec[1,2])

            # get optimal solutions
            v_opt   = u_vec[1,2]
            d_f_opt = u_vec[1,1]

            #v = v_opt*sqrt(1+c1*c1*d_f_opt*d_f_opt)
            w = -v_opt/l*d_f_opt
            car_control_msg = Twist2DStamped()
            car_control_msg.v = v_opt
            car_control_msg.omega = w
            #print(v)
            #publish commands
            publish(pub_car_cmd, car_control_msg)
        else
            println("unknown position")
        end
        rossleep(loop_rate)
    end
    spin()
end

if ! isinteractive()
    main()
end


text:
export poisson, poisson!

 # get the maximum intensity of an image
mymax(X::AbstractArray{<:RGB}) = max_rgb(X)
mymax(X::AbstractArray{<:Gray}) = gray(maximum(X))
mymax(X::AbstractArray) = maximum(X)

 # noise function for poisson
function noise_f(x, scaling, max_intens)
    # scale image to max_intensity and apply poisson noise
    # after Poisson noise scale it back
    return pois_rand(x * scaling / max_intens) * max_intens / scaling
end

f_pois(scaling, max_intens) = x -> scaling == nothing ? pois_rand(x) : noise_f(x, scaling, max_intens)

comb_pois(x, n) = n

poisson(X::AbstractArray, scaling=nothing; clip=false) = poisson!(copy(X), scaling, clip=clip)
function poisson!(X::AbstractArray, scaling=nothing; clip=false)
    max_intens = convert(Float64, mymax(X))
    return apply_noise!(comb_pois, f_pois(scaling, max_intens), X, clip)
end

"""
    poisson(X; scaling=nothing, clip=false)

Returns the array `X` affected by Poisson noise.
At every position the Poisson noise affects the intensity individually
and the values at the positions represent the expected value of the Poisson
Distribution.
Since Poisson Noise due to discrete events you should
provide the optional argument `scaling`. This `scaling` connects
the highest value of the array with the discrete number of events.
The highest value will be then scaled and the poisson noise is applied
Afterwards we scale the whole array back so that the initial intensity
is preserved but with applied Poisson noise.
`clip` is a keyword argument. If given, it clips the values to [0, 1]
"""
poisson

max_issues_repo_name: roflmaostc/Noise

max_forks_repo_name: roflmaostc/Noise


text:
using DecisionMakingProblems
using Test
using Random
using LinearAlgebra
using GridInterpolations

const p = DecisionMakingProblems

# MDP

@testset "2048.jl" begin
    m = TwentyFortyEight()
    mdp = MDP(m)
    @test length(mdp.𝒜) == 4
    @test mdp.γ == 1.0
    init_state = p.initial_board()
    s′, r = mdp.TR(init_state, p.DOWN)
    @test s′ != init_state || isdone(s′)
    @test r >= -1.0
end

@testset "cart_pole.jl" begin
    # m = CartPole(1.0, 10.0, 1.0, 1.0, 0.1, 9.8, 0.02, 4.8, deg2rad(24))
    m = CartPole()
    @test p.n_actions(m) == 2
    @test p.discount(m) == 1.0
    @test p.ordered_actions(m) == 1:2
    min_state = [-4.8, -Inf, -24 * pi/180, -Inf]
    max_state = [4.8, Inf, 24 * pi/180, Inf]
    state = p.generate_start_state(m)
    @test min_state <= p.vec(state) <= max_state
    @test !p.is_terminal(m, state)
    @test min_state <= p.vec(p.cart_pole_transition(m, state, rand(1:2))) <= max_state
    @test p.reward(m, state, rand(1:2)) in [0.0, 1.0]
    mdp = MDP(m)
end

@testset "collision_avoidance.jl" begin
    m = CollisionAvoidance()
    distrib = p.CollisionAvoidanceStateDistribution()
    s = p.rand(distrib)
    simple_pol = p.SimpleCollisionAvoidancePolicy()
    optimal_pol = p.OptimalCollisionAvoidancePolicy()
    @test length(p.vec(rand(p.transition(m, s, optimal_pol(s))))) == 4
    @test p.is_terminal(m, s) == (p.vec(s)[4] < 0.0)
    @test p.reward(m, rand(p.transition(m, s, optimal_pol(s))), rand(m.𝒜)) <= 0
    policy = p.CollisionAvoidanceValueFunction(m, simple_pol)
    mdp = MDP(m)
end

@testset "hexworld.jl" begin
    m = HexWorld()
    hexes = m.hexes
    @test p.n_states(m) == length(hexes) + 1 && p.ordered_states(m) == 1:length(hexes) + 1
    @test p.n_actions(m) == 6 && p.ordered_actions(m) == 1:6
    @test p.discount(m) == 0.9
    @test p.state_index(m, 1) == 1
    state = rand(1:p.n_states(m))
    action = rand(1:p.n_actions(m))
    p.transition(m, state, action)
    state_ = p.generate_s(m, state, action)
    reward = p.reward(m, state, action)
    @test state_ in p.ordered_states(m)
    @test reward <= 10
    @test p.generate_sr(m, state, action)[1] in p.ordered_states(m) && p.generate_sr(m, state, action)[2] <= 10
    @test p.generate_start_state(m) in p.ordered_states(m)
    @test p.hex_distance(rand(hexes), rand(hexes)) >= 0
    mdp = MDP(m)
end
@testset "simple_lqr.jl" begin
    m = LQR()
    @test p.discount(m) == 1.0
    state = p.generate_start_state(m)
    @test -10 <= rand(p.transition(m, state, rand())) <= 10
    @test p.reward(m, state, rand()) <= 0
    mdp = MDP(m)
end

@testset "mountain_car.jl" begin
    m = MountainCar()
    @test p.n_actions(m) == 3 && p.ordered_actions(m) == [1, 2, 3]
    @test p.discount(m) == 1.0
    state_min = [-1.2, -0.07]
    state_max = [0.6, 0.07]
    start_state = p.generate_start_state(m)
    @test all(state_min <= start_state <= state_max)
    @test all(state_min <= p.mountain_car_transition(start_state, 1) <= state_max)
    @test p.reward(m, start_state, 1) <= 0
    mdp = MDP(m)
end

# POMDP

@testset "crying_baby.jl" begin
    m = CryingBaby()
    @test p.n_states(m) == 2 && p.ordered_states(m) == [1, 2]
    @test p.n_actions(m) == 3 && p.ordered_actions(m) == [1, 2, 3]
    @test p.n_observations(m) == 2 && p.ordered_observations(m) == [true, false]
    @test p.discount(m) == 0.9
    p.transition(m, rand(1:2), rand(1:3))
    @test 0 <= p.observation(m, rand(1:3), rand(1:2)).p <= 1
    @test p.reward(m, rand(1:2), rand(1:3)) <= 0
    @test p.reward(m, [0.1, 0.9], rand(1:3)) <= 0
    pomdp = p.POMDP(m)
end

@testset "machine_replacement.jl" begin
    mdp = MachineReplacement()
    m = DiscretePOMDP(mdp)
    @test p.n_states(m) == 3 && p.ordered_states(m) == 1:3
    @test p.n_actions(m) == 4 && p.ordered_actions(m) == 1:4
    @test p.n_observations(m) == 2 && p.ordered_observations(m) == 1:2
    @test p.discount(m) == 1.0
    @test rand(p.transition(m, rand(1:3), rand(1:4))) in 1:3
    @test rand(p.observation(m, rand(1:4), rand(1:3))) in 1:2
    @test p.reward(m, rand(1:3), rand(1:4)) <= 1.0
    pomdp = p.POMDP(m)
end

@testset "catch.jl" begin
    mdp = Catch()
    m = DiscretePOMDP(mdp)
    @test p.n_states(m) == 4 && p.ordered_states(m) == 1:4
    @test p.n_actions(m) == 10 && p.ordered_actions(m) == 1:10
    @test p.n_observations(m) == 2 && p.ordered_observations(m) == 1:2
    @test p.discount(m) == 0.9
    @test rand(p.transition(m, rand(1:4), rand(1:10))) in 1:4
    @test rand(p.observation(m, rand(1:10), rand(1:4))) in 1:2
    @test p.reward(m, rand(1:4), rand(1:10)) >= 0
    pomdp = p.POMDP(m)
end

# Simple Game

@testset "prisoners_dilemma.jl" begin
    m = PrisonersDilemma()
    @test p.n_agents(m) == 2
    @test length(p.ordered_actions(m, rand(1:2))) == 2 && length(p.ordered_joint_actions(m)) == 4
    @test p.n_actions(m, rand(1:2)) == 2 && p.n_joint_actions(m) == 4
    @test p.reward(m, rand(1:2), [rand(p.ordered_actions(m, 0)), rand(p.ordered_actions(m, 0))]) <= 0.0
    @test p.joint_reward(m, [rand(p.ordered_actions(m, 0)), rand(p.ordered_actions(m, 0))]) <= [0.0, 0.0]
    simplegame = p.SimpleGame(m)
end

@testset "rock_paper_scissors.jl" begin
    m = RockPaperScissors()
    @test p.n_agents(m) == 2
    @test length(p.ordered_actions(m, rand(1:2))) == 3 && length(p.ordered_joint_actions(m)) == 9
    @test p.n_actions(m, rand(1:2)) == 3 && p.n_joint_actions(m) == 9
    @test -1.0 <= p.reward(m, rand(1:2), [rand(p.ordered_actions(m, 0)), rand(p.ordered_actions(m, 0))]) <= 1.0
    @test [-1.0, -1.0] <= p.joint_reward(m, [rand(p.ordered_actions(m, 0)), rand(p.ordered_actions(m, 0))]) <= [1.0, 1.0]
    simplegame = p.SimpleGame(m)
end

@testset "travelers.jl" begin
    m = Travelers()
    @test p.n_agents(m) == 2
    @test length(p.ordered_actions(m, rand(1:2))) == 99 && length(p.ordered_joint_actions(m)) == 99^2
    @test p.n_actions(m, rand(1:2)) == 99 && p.n_joint_actions(m) == 99^2
    @test 0.0 <= p.reward(m, rand(1:2), [rand(p.ordered_actions(m, 0)), rand(p.ordered_actions(m, 0))]) <= 102
    @test [0.0, 0.0] <= p.joint_reward(m, [rand(p.ordered_actions(m, 0)), rand(p.ordered_actions(m, 0))]) <= [102, 102]
    simplegame = p.SimpleGame(m)
end

# Markov Game

@testset "predator_prey.jl" begin
    m = PredatorPreyHexWorld()
    hexes = m.hexes
    @test p.n_agents(m) == 2
    @test length(p.ordered_states(m, rand(1:2))) == length(hexes) && length(p.ordered_states(m)) == length(hexes)^2
    @test length(p.ordered_actions(m, rand(1:2))) == 6 && length(p.ordered_joint_actions(m)) == 36
    @test p.n_actions(m, rand(1:2)) == 6 && p.n_joint_actions(m) == 36

    @test 0.0 <= p.transition(m, rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m)), rand(p.ordered_states(m))) <= 1.0
    @test -1.0 <= p.reward(m, rand(1:2), rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m))) <= 10.0
    @test [-1.0, -1.0] <= p.joint_reward(m, rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m))) <= [10.0, 10.0]
    mg = p.MG(m)
end

# POMG

@testset "multicaregiver.jl" begin
    m = MultiCaregiverCryingBaby()
    @test p.n_agents(m) == 2
    @test length(p.ordered_states(m)) == 2
    @test length(p.ordered_actions(m, rand(1:2))) == 3 && length(p.ordered_joint_actions(m)) == 9
    @test p.n_actions(m, rand(1:2)) == 3 && p.n_joint_actions(m) == 9
    @test length(p.ordered_observations(m, rand(1:2))) == 2 && length(p.ordered_joint_observations(m)) == 4
    @test p.n_observations(m, rand(1:2)) == 2 && p.n_joint_observations(m) == 4

    @test 0.0 <= p.transition(m, rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m)), rand(p.ordered_states(m))) <= 1.0
    @test 0.0 <= p.joint_observation(m, rand(p.ordered_joint_actions(m)), rand(p.ordered_states(m)), rand(p.ordered_joint_observations(m))) <= 1.0
    @test p.joint_reward(m, rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m))) <= [0.0, 0.0]
    @test p.joint_reward(m, rand(Float64, 2), rand(p.ordered_joint_actions(m))) <= [0.0, 0.0]
    pomg = p.POMG(m)
end

# DecPOMDP

@testset "collab_predator_prey.jl" begin
    m = CollaborativePredatorPreyHexWorld()
    hexes = m.hexes
    @test p.n_agents(m) == 2
    @test length(p.ordered_states(m, rand(1:2))) == length(hexes) && length(p.ordered_states(m)) == length(hexes)^3
    @test length(p.ordered_actions(m, rand(1:2))) == 6 && length(p.ordered_joint_actions(m)) == 36
    @test p.n_actions(m, rand(1:2)) == 6 && p.n_joint_actions(m) == 36
    @test length(p.ordered_observations(m, rand(1:2))) == 2 && length(p.ordered_joint_observations(m)) == 4
    @test p.n_observations(m, rand(1:2)) == 2 && p.n_joint_observations(m) == 4

    @test 0.0 <= p.transition(m, rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m)), rand(p.ordered_states(m))) <= 1.0
    @test 0.0 <= p.joint_observation(m, rand(p.ordered_joint_actions(m)), rand(p.ordered_states(m)), rand(p.ordered_joint_observations(m))) <= 1.0
    @test -1.0 <= p.reward(m, rand(p.ordered_states(m)), rand(p.ordered_joint_actions(m))) <= 10.0
    decpomdp = p.DecPOMDP(m)
end


text:
function solve()
    l = parse(Int, readline())
    ret = 1
    for (x,y)=zip(l-11:l-1, 1:11)
        ret *= x
        ret ÷= y
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# Testing the function
using Test
using Feather

using RegressionDiscontinuity
using Distributions
using GLM
using DataFrames

# Testing that the package doesn't detect false positives:
R1 = rand(Normal(), 100_000)
rv1 = RunningVariable(R1)
test1 = fit(McCraryTest(), rv1)

R2 = rand(Uniform(-1, 1), 100_000)
rv2 = RunningVariable(R2)
test2 = fit(McCraryTest(), rv2)

#@test test1.pval >= 0.05
#@test test2.pval >= 0.05

# Testing that the package detect true positives
R3 = rand(Uniform(-1, 1), 100_000)
R3 = R3 .+ 2 * (1.0 * (rand(Uniform(-1, 1), 100_000) .> 0) .* (R3 .< 0))
rv3 = RunningVariable(R3)
test3 = fit(McCraryTest(), rv3)

#@test test3.pval < 0.05

# Testing the package with real data.
lee08 = load_rdd_data(:lee08)
lee08_mccrary = fit(McCraryTest(), lee08.ZsR)

# I will compare the results with the rdd package.
pvalR = 0.1982771
binR = 0.011
bwR = 0.242
θhatR = 0.103

@test lee08_mccrary.pval ≈ pvalR atol = 0.001
@test lee08_mccrary.b ≈ binR atol = 0.001
@test lee08_mccrary.h ≈ bwR atol = 0.001
@test lee08_mccrary.θhat ≈ θhatR atol = 0.001


text:
import QuadratureRules: GaussLegendreQuadrature, LobattoLegendreQuadrature

"""
Lagrange basis on the interval [0..1].
"""
struct Lagrange{T, BT, XT <: AbstractVector{T}} <: Basis{T}
    b::BT
    x::XT

    denom::XT
    diffs::Matrix{T}
    vdminv::Matrix{T}

    function Lagrange{T}(x::XT) where {T, XT <: SVector}
        n = length(x)
        denom = zeros(n)
        diffs = zeros(n,n)

        for i in eachindex(x)
            local p = 1
            for j in eachindex(x)
                diffs[i,j] = x[i] - x[j]
                if i ≠ j
                    p *= diffs[i,j]
                end
            end
            denom[i] = 1/p
        end

        sdenom = SVector{n}(denom)

        b = collect(y -> sdenom[j] * mapreduce(i -> i ≠ j ? (y - x[i]) : one(T), *, eachindex(x)) for j in eachindex(sdenom))

        new{T, typeof(b), typeof(x)}(b, x, sdenom, diffs, vandermonde_matrix_inverse(x))
    end

    Lagrange{T}(x::Vector) where {T} = Lagrange{T}(SVector{length(x),T}(x))
    Lagrange{T}(x::AbstractVector) where {T} = Lagrange{T}(collect(x))
end

Lagrange(x::AbstractVector{T}) where {T} = Lagrange{T}(x)

LagrangeGauß(n) = Lagrange(GaussLegendreQuadrature(n).nodes)
LagrangeLobatto(n) = Lagrange(LobattoLegendreQuadrature(n).nodes)

(L::Lagrange)(x::Number, j::Integer) = L.b[j](x)

basis(L::Lagrange) = L.b
nodes(L::Lagrange) = L.x
nbasis(L::Lagrange) = length(basis(L))
nnodes(L::Lagrange) = length(nodes(L))
order(L::Lagrange)  = nnodes(L)
degree(L::Lagrange) = nnodes(L) - 1

Base.eltype(::Lagrange{T}) where {T} = T
Base.eachindex(L::Lagrange) = eachindex(L.b)
Base.axes(L::Lagrange) = (Inclusion(0..1), eachindex(L))
ContinuumArrays.grid(L::Lagrange) = nodes(L)

Base.hash(L::Lagrange, h::UInt) = hash(L.x, h)
Base.:(==)(L1::Lagrange, L2::Lagrange) = (L1.x == L2.x)
Base.isequal(L1::Lagrange{T1}, L2::Lagrange{T2}) where {T1,T2} = (T1 == T2 && L1 == L2)
Base.isapprox(L1::Lagrange, L2::Lagrange; kwargs...) = isapprox(L1.x, L2.x; kwargs...)

Base.getindex(L::Lagrange, x::Number, j::Integer) = L(x,j)
Base.getindex(L::Lagrange, x::Number,  ::Colon) = [b(x) for b in L.b]
Base.getindex(L::Lagrange, X::AbstractVector, j::Integer) = L.(X,j)
Base.getindex(L::Lagrange, X::AbstractVector,  ::Colon) = [b(x) for x in X, b in L.b]

## Derivative

@simplify *(D::Derivative, L::Lagrange) = Mul(D,L)

const LagrangeDerivative = QMul2{<:Derivative,<:Lagrange}

function _eval(D::LagrangeDerivative, x::DT, j::Int) where {DT}
    local L = D.B
    local T = promote_type(eltype(L), DT)
    local d::T = 0

    for l in eachindex(L)
        if l ≠ j
            z = 1 / L.diffs[j,l]
            for i in eachindex(L)
                i ≠ j && i ≠ l ? z *= (x - L.x[i]) / L.diffs[j,i] : nothing
            end
            d += z
        end
    end
    return d
end

Base.getindex(D::LagrangeDerivative, x::Number, j::Integer) = _eval(D, x, j)
Base.getindex(D::LagrangeDerivative, x::Number,  ::Colon) = [_eval(D, x, j) for j in eachindex(D.B)]
Base.getindex(D::LagrangeDerivative, X::AbstractVector, j::Integer) = [_eval(D, x, j) for x in X]
Base.getindex(D::LagrangeDerivative, X::AbstractVector,  ::Colon) = [_eval(D, x, j) for x in X, j in eachindex(D.B)]

Base.adjoint(L::Lagrange) = Derivative(axes(L,1)) * L


text:
using JSON

JSON_PATH = joinpath(dirname(@__FILE__), "../data/hr/hr.json")
citydb = JSON.parsefile(JSON_PATH)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
function send_command(scope, cmd, args...)
    scope["_commands"][] = [cmd, args...]
    # The handler for _commands is set up when plot is constructed
    nothing
end

# ----------------------- #
# Plotly.js api functions #
# ----------------------- #

# for each of these we first update the Julia object, then update the display

function restyle!(
        plt::WebIOPlot, ind::Union{Int,AbstractVector{Int}},
        update::AbstractDict=Dict();
        kwargs...)
        restyle!(plt.p, ind, update; kwargs...)
        send_command(plt.scope, :restyle, merge(update, prep_kwargs(kwargs)), ind .- 1)
end

function restyle!(plt::WebIOPlot, update::AbstractDict=Dict(); kwargs...)
    restyle!(plt.p, update; kwargs...)
    send_command(plt.scope, :restyle, merge(update, prep_kwargs(kwargs)))
end

function relayout!(plt::WebIOPlot, update::AbstractDict=Dict(); kwargs...)
    relayout!(plt.p, update; kwargs...)
    send_command(plt.scope, :relayout, merge(update, prep_kwargs(kwargs)))
end

function update!(
        plt::WebIOPlot, ind::Union{Int,AbstractVector{Int}},
        update::AbstractDict=Dict();
        layout::Layout=Layout(),
        kwargs...)
    update!(plt.p, ind, update; layout=layout, kwargs...)
    send_command(plt.scope, :update, merge(update, prep_kwargs(kwargs)), layout, ind .- 1)
end

function update!(
        plt::WebIOPlot, update::AbstractDict=Dict(); layout::Layout=Layout(),
        kwargs...
    )
    update!(plt.p, update; layout=layout, kwargs...)
    send_command(plt.scope, :update, merge(update, prep_kwargs(kwargs)), layout)
end

function addtraces!(plt::WebIOPlot, traces::AbstractTrace...)
    addtraces!(plt.p, traces...)
    send_command(plt.scope, :addTraces, traces)
end

function addtraces!(plt::WebIOPlot, i::Int, traces::AbstractTrace...)
    addtraces!(plt.p, i, traces...)
    send_command(plt.scope, :addTraces, traces, i-1)
end

function deletetraces!(plt::WebIOPlot, inds::Int...)
    deletetraces!(plt.p, inds...)
    send_command(plt.scope, :deleteTraces, collect(inds) .- 1)
end

function movetraces!(plt::WebIOPlot, to_end::Int...)
    movetraces!(plt.p, to_end...)
    send_command(plt.scope, :moveTraces, traces, collect(to_end) .- 1)
end

function movetraces!(
        plt::WebIOPlot, src::AbstractVector{Int}, dest::AbstractVector{Int}
    )
    movetraces!(plt.p, src, dest)
    send_command(plt.scope, :moveTraces, traces, src .- 1, dest .- 1)
end

function redraw!(plt::WebIOPlot)
    redraw!(plt.p)
    send_command(plt.scope, :redraw)
end

function purge!(plt::WebIOPlot)
    purge!(plt.p)
    send_command(plt.scope, :purge)
end

function to_image(plt::WebIOPlot; kwargs...)
    to_image(plt.p)
    send_command(plt.scope, :toImage, Dict(kwargs))
end

function download_image(plt::WebIOPlot; kwargs...)
    download_image(plt.p)
    send_command(plt.scope, :downloadImage, Dict(kwargs))
end

# unexported (by plotly.js) api methods
function extendtraces!(plt::WebIOPlot, update::AbstractDict=Dict(),
              indices::AbstractVector{Int}=[1], maxpoints=-1;)
    extendtraces!(plt.p, update, indices, maxpoints)
    send_command(
        plt.scope, :extendTraces, prep_kwargs(update), indices .- 1, maxpoints
    )
end

function prependtraces!(plt::WebIOPlot, update::AbstractDict=Dict(),
               indices::AbstractVector{Int}=[1], maxpoints=-1;)
    prependtraces!(plt.p, update, indices, maxpoints)
    send_command(
        plt.scope, :prependTraces, prep_kwargs(update), indices .- 1, maxpoints
    )
end

for f in [:restyle, :relayout, :update, :addtraces, :deletetraces,
          :movetraces, :redraw, :extendtraces, :prependtraces, :purge]
    f! = Symbol(f, "!")
    @eval function $(f)(plt::WebIOPlot, args...; kwargs...)
        out = WebIOPlot(deepcopy(plt.p))
        $(f!)(out, args...; kwargs...)
        out
    end
end


text:
@inline zerofunction(args...) = 0
@inline onefunction(args...) = 1

"""
    struct RelaxingFunction{R, M, T}

Callable object for restoring fields to a `target` at
some `rate` and within a `mask`ed region in `x, y, z`.
"""
struct RelaxingFunction{R, M, T}
      rate :: R
      mask :: M
    target :: T
end

@inline (f::RelaxingFunction)(x, y, z, t, field) =
    f.rate * f.mask(x, y, z) * (f.target(x, y, z, t) - field)

#####
##### Relaxing forcing
#####

"""
    Relaxation(; rate, mask=onefunction, target=zerofunction)

Returns a `SimpleForcing` that restores a field to `target(x, y, z, t)`
at the specified `rate`, in the region `mask(x, y, z)`.

Example
=======

* Restore a field to a linear z-gradient everywhere on a timescale of "60" (equal
  to one minute if the time units of the simulation are seconds).

```julia
julia> restore_stratification = Relaxation(; rate = 1/60, target = (x, y, z, t) -> z)
```
"""
Relaxation(; rate, mask=onefunction, target=zerofunction) =
    SimpleForcing(RelaxingFunction(rate, mask, target); field_in_signature=true)

#####
##### Sponge layer functions
#####

"""
    GaussianMask{D}(center, width)

Callable object that returns a Gaussian masking function centered on
`center`, with `width`, and varying along direction `D`.

Examples
========

* Create a Gaussian mask centered on `z=0` with width `1` meter.

```julia
julia> mask = GaussianMask{:z}(0, 1)
```
"""
struct GaussianMask{D, T}
    center :: T
     width :: T

    function GaussianMask{D}(; center, width) where D
        T = promote_type(typeof(center), typeof(width))
        return new{D, T}(center, width)
    end
end

@inline (g::GaussianMask{:x})(x, y, z) = exp(-(x - g.center)^2 / (2 * g.width^2))
@inline (g::GaussianMask{:y})(x, y, z) = exp(-(y - g.center)^2 / (2 * g.width^2))
@inline (g::GaussianMask{:z})(x, y, z) = exp(-(z - g.center)^2 / (2 * g.width^2))

#####
##### Linear target functions
#####

"""
    LinearTarget{D}(intercept, gradient)

Callable object that returns a Linear target function
with `intercept` and `gradient`, and varying along direction `D`.

Examples
========

* Create a linear target function varying in `z`, equal to `0` at
  `z=0` and with gradient 10⁻⁶:

```julia
julia> target = LinearTarget{:z}(0, 1e-6)
```
"""
struct LinearTarget{D, T}
    intercept :: T
     gradient :: T

    function LinearTarget{D}(; intercept, gradient) where D
        T = promote_type(typeof(gradient), typeof(intercept))
        return new{D, T}(intercept, gradient)
    end
end

@inline (p::LinearTarget{:x})(x, y, z, t) = p.intercept + p.gradient * x
@inline (p::LinearTarget{:y})(x, y, z, t) = p.intercept + p.gradient * y
@inline (p::LinearTarget{:z})(x, y, z, t) = p.intercept + p.gradient * z