text:
struct SimpleUpdateProblem{TU,TA,TW}
    u::TU
    a::TA
    b::TA
    σab::TW
    σba::TW
    ωab::TW
    ωba::TW
end

"""
    simpleupdate(u, ainit, binit)
creates an object of type `SimpleUpdateProblem` that contains the propagator `u`,
the intial tensors `a` and `b` as well as the weights σab, σba, ωab and ωba.
Graphically the structure looks like:

          |_________|
          [____u____]
          |         |
          |  ωab    | ωba
          | /       | /
    -σba-[a]--σab--[b]-σba-
         /         /
       ωba       ωab

with the indices-ordering like:

    3|________4|
     [____u____]
    1|        2|

       5| /3
    4--[a]--2
       /1

"""
function simpleupdate(u, ainit, binit)
    a, b = deepcopy(ainit), deepcopy(binit)
    σab = TensorOperations.similar_from_indices(
        real(eltype(a)), (2,), (4,), (1,),(2,), a, b, :Y, :Y)
    σba, ωab, ωba = deepcopy(σab), deepcopy(σab), deepcopy(σab)
    apply!.((σab, σba, ωab, ωba), x -> x .= Matrix(I, size(x)...))
    return SimpleUpdateProblem{typeof(u), typeof(a), typeof(σab)}(u, a, b, σab, σba, ωab, ωba)
end

update(su::SimpleUpdateProblem) = update!(deepcopy(su))

function update!(su::SimpleUpdateProblem)
    @unpack a, b, σab, σba, ωab, ωba, u = su

    update!((a,b,σab), (σba, ωab, ωba, u))
    update!((b,a,σba), (σab, ωba, ωab, u))

    @tensor begin
        arot[1,2,3,4,5] := a[4,1,2,3,5]
        brot[1,2,3,4,5] := b[4,1,2,3,5]
    end

    update!((arot, brot, ωab), (ωba, σab, σba, u))
    update!((brot, arot, ωba), (ωab, σba, σab, u))

    @tensor begin
        a[1,2,3,4,5] = arot[2,3,4,1,5]
        b[1,2,3,4,5] = brot[2,3,4,1,5]
    end

    return su
end

function update!((a,b,σab), (σba, ωab, ωba, u))
    #contract with weights
    @tensor begin
        aw[1,2,3,4,5] := σba[4,-4] * a[-1,2,-3,-4,5] * ωba[1,-1] * ωab[-3,3]
        bw[1,2,3,4,5] := σba[-2,2] * b[-1,-2,-3,4,5] * ωab[1,-1] * ωba[-3,3]
    end

    #reduce tensors and apply u
    qa, ra = tensorqr(aw, ((1,3,4),(5,2)))
    rb, qb = tensorrq(bw, ((4,5),(1,2,3)))

    @tensor θ[1,2,3,4] := ra[1,p1,-1] * σab[-1,-2] * rb[-2,p2,3] * u[p1,p2,2,4]

    #truncate bond dim
    s = a isa DTensor ? size(a,2) : sizes(a,2)[1]
    Uθ, Σθ, Vdθ = tensorsvd(θ, ((1,2),(3,4)), svdtrunc = svdtrunc_maxχ(s))

    @tensor begin
        aw[1,2,3,4,5] = qa[1,3,4,-1] * Uθ[-1,5,2]
        bw[1,2,3,4,5] = Vdθ[4,-4,5] * qb[-4,1,2,3]
    end

    #reextract weights
    iσab, iσba, iωab, iωba = pinv.((σab, σba, ωab, ωba))

    @tensor begin
        a[1,2,3,4,5] = iσba[4,-4] * aw[-1,2,-3,-4,5] * iωba[1,-1] * iωab[-3,3]
        b[1,2,3,4,5] = iσba[-2,2] * bw[-1,-2,-3,4,5] * iωab[1,-1] * iωba[-3,3]
    end

     apply!(Σθ, x -> x .= x ./ sum(diag(x)))
     copyto!(σab, Σθ)
     return a, b, σab
end

"""
    energy(su::SimpleUpdateProblem, h)
returns the energy calculated with diagrams such as:

            √ωab      √ωba
            / |       /|
    √σba-[a*]-σab--[b*]√σba
    |    /|   |    /|  |  |
    | √ωba|   | √ωab|  |  |
    |   | |___|___|_|  |  |
    |   | [____h__|_]  |  |
    |   | | √ωab  | |√ωba |
    |   | | /     | | /   |
    √σba|[a]--σab-|[b]√σba-
        |/        |/
      √ωba      √ωab

which might not be the correct approach
"""
function energy(su::SimpleUpdateProblem, h)
    @unpack a, b, σab, σba, ωab, ωba = su
    @tensor aw[1,2,3,4,5] := a[-1,-2,-3,-4,5] *
                             ωba[1,-1] *
                             σab[-2,2] *
                             ωab[-3,3] *
                             σba[4,-4]

    #E horizontal a-b
    @tensor bw[1,2,3,4,5] := b[-1,-2,-3,4,5] *
                             ωab[1,-1] *
                             σba[2,-2] *
                             ωba[-3,3]
    @tensor ehab = aw[-1,l1,-3,-4,1] * conj(aw)[-1,l2,-3,-4,3] *
                   bw[-5,-6,-7,l1,2] * conj(bw)[-5,-6,-7,l2,4] *
                   h[1,2,3,4]
    @tensor nhab = aw[-1,l1,-3,-4,i1] * conj(aw)[-1,l2,-3,-4,i1] *
                   bw[-5,-6,-7,l1,i2] * conj(bw)[-5,-6,-7,l2,i2]

    #E horizontal b-a
    @tensor bw[1,2,3,4,5] := b[-1,2,-3,-4,5] *
                             ωab[1,-1] *
                             ωba[-3,3] *
                             σab[4,-4]
    @tensor ehba = bw[-1,l1,-3,-4,1] * conj(bw)[-1,l2,-3,-4,3] *
                   aw[-5,-6,-7,l1,2] * conj(aw)[-5,-6,-7,l2,4] *
                   h[1,2,3,4]
    @tensor nhba = bw[-1,l1,-3,-4,i1] * conj(bw)[-1,l2,-3,-4,i1] *
                   aw[-5,-6,-7,l1,i2] * conj(aw)[-5,-6,-7,l2,i2]

    #E vertical a-b
    @tensor bw[1,2,3,4,5] := b[1,-2,-3,-4,5] *
                             σba[2,-2] *
                             ωba[-3,3] *
                             σab[4,-4]
    @tensor evab = aw[-1,-2,l1,-4,1] * conj(aw)[-1,-2,l2,-4,3] *
                   bw[l1,-6,-7,-8,2] * conj(bw)[l2,-6,-7,-8,4] *
                   h[1,2,3,4]
    @tensor nvab = aw[-1,-2,l1,-4,i1] * conj(aw)[-1,-2,l2,-4,i1] *
                   bw[l1,-6,-7,-8,i2] * conj(bw)[l2,-6,-7,-8,i2]

    #E vertical b-a
    @tensor bw[1,2,3,4,5] := b[-1,-2,3,-4,5] *
                             ωab[1,-1] *
                             σba[2,-2] *
                             σab[4,-4]
    @tensor evba = bw[-1,-2,l1,-4,1] * conj(bw)[-1,-2,l2,-4,3] *
                   aw[l1,-6,-7,-8,2] * conj(aw)[l2,-6,-7,-8,4] *
                   h[1,2,3,4]
    @tensor nvba = bw[-1,-2,l1,-4,i1] * conj(bw)[-1,-2,l2,-4,i1] *
                   aw[l1,-6,-7,-8,i2] * conj(aw)[l2,-6,-7,-8,i2]
    es = (ehab/nhab,  ehba/nhba, evab/nvab, evba/nvba)
    return Real.(es)
    return sum(es)/2
end


text:
function pisum()
	sum = 0.0
	for j = 1:500
		sum = 0.0
		for k = 1:10000
			sum += 1.0/(k*k)
		end
	end
	sum
end

start = time()
temp = pisum()
stop = time()
elapsed = (stop - start)

println(elapsed)

ext: jl

max_stars_repo_name: owainkenwayucl/Garbage

max_issues_repo_name: owainkenwayucl/Garbage

max_forks_repo_name: owainkenwayucl/Garbage


text:
# TO DO: [200], [400], [405]
# [200] Generic Event Detection Blockette (52 bytes)
# [400] Beam Blockette (16 bytes)
# [405] Beam Delay Blockette (6 bytes)

# All blockette functions return a UInt16 value equal to blockette length in bytes
function blk_time!(t::Array{Int32,1}, sid::IO, b::Bool)
  yy    = read(sid, UInt16)
  jj    = read(sid, UInt16)
  t[4]  = Int32(read(sid, UInt8))
  t[5]  = Int32(read(sid, UInt8))
  t[6]  = Int32(read(sid, UInt8))
  skip(sid, 1)
  ms    = read(sid, UInt16)
  if b
    yy = ntoh(yy)
    jj = ntoh(jj)
    ms = ntoh(ms)
  end
  yy = Int32(yy)
  jj = Int32(jj)
  t[1] = yy
  (t[2], t[3]) = j2md(yy, jj)
  t[7] = div(Int32(ms), Int32(10))
  return nothing
end

# [100] Sample Rate Blockette (12 bytes)
function blk_100(S::SeisData, sid::IO, c::Int64)
  BUF.dt = Float64(BUF.swap ? ntoh(read(sid, Float32)) : read(sid, Float32))
  skip(sid, 4)
  return 0x000c
end

# [201] Murdock Event Detection Blockette (60 bytes)
function blk_201(S::SeisData, sid::IO, c::Int64)
  read!(sid, BUF.B201.sig)                              # amplitude, period, background estimate
  BUF.B201.flags = read(sid, UInt8)                     # detection flags
  skip(sid, 1)                                          # reserved
  blk_time!(BUF.B201.t, sid, BUF.swap)                  # onset time
  read!(sid, BUF.B201.snr)                              # snr, lookback, pick algorithm
  BUF.B201.det  = String(read(sid, 24))                 # detector name

  # will store as a string: t_evt, :sig (3 vals), bittsring(:flags), :snr, :lb, :p)
  # Store event detections in S.misc
  sig = BUF.B201.sig
  flag = BUF.B201.flags == 0x80 ? "dilatation" : "compression"
  if BUF.swap
    sig = ntoh.(sig)
    flag = BUF.B201.flags == 0x01 ? "dilatation" : "compression"
  end
  if !haskey(S.misc[c], "seed_event")
    S.misc[c]["seed_event"] = Array{String, 1}(undef,0)
  end
  push!(S.misc[c]["seed_event"],  join(BUF.B201.t, ',') * "," *
                                    join(sig, ',') * "," *
                                    flag * "," *
                                    join(BUF.B201.snr, ',') * "," *
                                    strip(BUF.B201.det) )
  return 0x003c
end

# [300] Step Calibration Blockette (60 bytes)
# [310] Sine Calibration Blockette (60 bytes)
# [320] Pseudo-random Calibration Blockette (64 bytes)
# [390] Generic Calibration Blockette (28 bytes)
function blk_calib(S::SeisData, sid::IO, c::Int64, bt::UInt16)
  p = position(sid)
  blk_time!(BUF.Calib.t, sid, BUF.swap)       # Calibration time
  skip(sid, 1)                                  # Reserved byte
  if bt == 0x012c
    BUF.Calib.n        = read(sid, UInt8)      # Number of step calibrations
  end
  BUF.Calib.flags      = read(sid, UInt8)      # Calibration flags
  BUF.Calib.dur1       = read(sid, UInt32)     # Calibration duration
  if bt == 0x012c
    BUF.Calib.dur2     = read(sid, UInt32)     # Interval duration
  elseif bt == 0x0136
    BUF.Calib.period   = read(sid, Float32)    # Period of signal (seconds)
  end
  BUF.Calib.amplitude  = read(sid, Float32)    # Peak-to-peak amplitude
  BUF.Calib.channel    = read(sid, 3)
  skip(sid, 1)                                  # Reserved byte
  BUF.Calib.ref        = read(sid, UInt32)     # Reference amplitude

  # String arrays
  if bt < 0x0186
    BUF.Calib.coupling = read(sid, 12)
    BUF.Calib.rolloff  = read(sid, 12)
    if bt == 0x0140
      BUF.Calib.noise  = read(sid, 8)
    end
  else
    BUF.Calib.coupling = Array{UInt8,1}(undef,0)
    BUF.Calib.rolloff  = Array{UInt8,1}(undef,0)
  end
  bc = Char[' ', '\0']

  # Check that we already have "seed_calib"
  if !haskey(S.misc[c], "seed_calib")
    S.misc[c]["seed_calib"] = Array{String, 1}(undef,0)
  end

  # Swap as needed
  if BUF.swap
    flag = reverse(bitstring(BUF.Calib.flags))
    amp = ntoh(BUF.Calib.amplitude)
    ref = ntoh(BUF.Calib.ref)
    per = ntoh(BUF.Calib.period)
    BUF.Calib.dur1 = ntoh(BUF.Calib.dur1)
    BUF.Calib.dur2 = ntoh(BUF.Calib.dur2)
  else
    flag = bitstring(BUF.Calib.flags)
    amp = BUF.Calib.amplitude
    ref = BUF.Calib.ref
    per = BUF.Calib.period
  end
  typ = (bt == 0x012c ? "Step" : (bt == 310 ? "Sine" : (bt == 320 ? "Pseudo-random" : "Generic")))
  dur = bt == 0x012c ? string(BUF.Calib.dur1 * "," * BUF.Calib.dur2) : string(BUF.Calib.dur1)
  amp = bt == 0x0136 ? join([string(per), string(amp)], ",") : string(amp)
  flag = bt == 0x012c ? flag * "," * string(BUF.Calib.n) : flag
  calib_str = join( [ typ,
                      join(BUF.Calib.t, ","),
                      flag,
                      dur,
                      amp,
                      strip(String(BUF.Calib.channel), bc),
                      string(ref) ], "," )
  if bt < 0x0186
    calib_str *= string(",", strip(String(BUF.Calib.coupling), bc), ",",
                        strip(String(BUF.Calib.rolloff), bc))
    if bt == 0x0140
      calib_str *= string(",", strip(String(BUF.Calib.noise), bc))
    end
  end
  push!(S.misc[c]["seed_calib"], calib_str)
  return UInt16(position(sid)-p) + 0x0004
end

# [500] Timing Blockette (200 bytes)
function blk_500(S::SeisData, sid::IO, c::Int64)
  BUF.B500.vco_correction    = BUF.swap ? ntoh(read(sid, Float32)) : read(sid, Float32)
  blk_time!(BUF.B500.t, sid, BUF.swap)
  BUF.B500.μsec              = read(sid, Int8)
  BUF.B500.reception_quality = read(sid, UInt8)
  BUF.B500.exception_count   = BUF.swap ? ntoh(read(sid, UInt32)) : read(sid, UInt32)
  BUF.B500.exception_type    = strip(String(read(sid, 16)), ['\0',' '])
  BUF.B500.clock_model       = strip(String(read(sid, 32)), ['\0',' '])
  BUF.B500.clock_status      = strip(String(read(sid, 128)), ['\0',' '])

  # TO DO: something with this
  if !haskey(S.misc[c], "seed_timing")
    S.misc[c]["seed_timing"] = Array{String, 1}(undef,0)
  end
  push!(S.misc[c]["seed_timing"], join([string(BUF.B500.vco_correction),
                                        join(BUF.B500.t, ','),
                                        BUF.B500.μsec,
                                        BUF.B500.reception_quality,
                                        BUF.B500.exception_count,
                                        BUF.B500.exception_type,
                                        BUF.B500.clock_model,
                                        BUF.B500.clock_status], ',')
        )
  return 0x00c8
end
# TO DO: correct S.t[c] when one of these timing blockettes is detected

# [1000] Data Only SEED Blockette (8 bytes)
function blk_1000(S::SeisData, sid::IO, c::Int64)
  BUF.fmt  = read(sid, UInt8)
  BUF.wo   = read(sid, UInt8)
  lx        = read(sid, UInt8)
  skip(sid, 1)
  BUF.nx   = 2^lx
  BUF.xs   = ((BUF.swap == true) && (BUF.wo == 0x01))
  return 0x0008
end

# [1001] Data Extension Blockette (8 bytes)
function blk_1001(S::SeisData, sid::IO, c::Int64)
  skip(sid, 1)
  BUF.tc += read(sid, Int8)
  skip(sid, 2)
  return 0x0008
end

# [2000] Variable Length Opaque Data Blockette
function blk_2000(S::SeisData, sid::IO, c::Int64)
  # Always big-Endian? Undocumented
  BUF.B2000.NB     = BUF.swap ? ntoh(read(sid, UInt16)) : read(sid, UInt16)
  BUF.B2000.os     = BUF.swap ? ntoh(read(sid, UInt16)) : read(sid, UInt16)
  n                 = BUF.swap ? ntoh(read(sid, UInt32)) : read(sid, UInt32)
  wo                = read(sid, UInt8)
  nf                = read(sid, UInt8)
  flag              = read(sid, UInt8)
  BUF.B2000.hdr    = read(sid, Int(BUF.B2000.os)-15)
  BUF.B2000.data   = read(sid, BUF.B2000.NB-BUF.B2000.os)

  # Store to S.misc[i]
  r = "seed_opaque_" * string(n)
  S.misc[c][r * "_wo"]    = wo
  S.misc[c][r * "_flag"]  = bitstring(BUF.B2000.flag)
  S.misc[c][r * "_hdr"]   = String.(split(String(BUF.B2000.hdr), '~', keepempty=true, limit=Int(nf)))
  S.misc[c][r * "_data"]  = BUF.B2000.data
  return BUF.B2000.NB
end


text: println("Use Shell to run.")


text:
function move_contributions(
  scell_to_val::AbstractArray, strian::Triangulation)

  scell_to_val, strian
end

function move_contributions(
  scell_to_val::AbstractArray, strian::Triangulation, ttrian::Triangulation)

  D = num_cell_dims(ttrian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  move_contributions(scell_to_val,sglue,tglue)
end

function move_contributions(
  scell_to_val::AbstractArray, sglue::FaceToFaceGlue, tglue::FaceToFaceGlue)
  move_contributions(scell_to_val,sglue,tglue,sglue.mface_to_tface)
end

function move_contributions(
  scell_to_val::AbstractArray,
  sglue::FaceToFaceGlue,
  tglue::FaceToFaceGlue,
  mcell_to_scell::AbstractArray)

  mcell_to_val = extend(scell_to_val,mcell_to_scell)
  tcell_to_mcell = tglue.tface_to_mface
  tcell_to_val = lazy_map(Reindex(mcell_to_val),tcell_to_mcell)
end

function move_contributions(
  scell_to_val::AbstractArray,
  sglue::FaceToFaceGlue,
  tglue::FaceToFaceGlue,
  mcell_to_scell::Nothing)

  scell_to_mcell = sglue.tface_to_mface
  mcell_to_tcell = tglue.mface_to_tface
  ntcells = length(tglue.tface_to_mface)
  tcell_to_scells_ptr = zeros(Int32,ntcells+1)
  for mcell in scell_to_mcell
    tcell = mcell_to_tcell[mcell]
    tcell_to_scells_ptr[1+tcell] += Int32(1)
  end
  length_to_ptrs!(tcell_to_scells_ptr)
  ndata = tcell_to_scells_ptr[end]-1
  tcell_to_scells_data = zeros(Int32,ndata)
  for (scell,mcell) in enumerate(scell_to_mcell)
    tcell = mcell_to_tcell[mcell]
    p = tcell_to_scells_ptr[tcell]
    tcell_to_scells_data[p] = scell
    tcell_to_scells_ptr[tcell] += Int32(1)
  end
  rewind_ptrs!(tcell_to_scells_ptr)
  tcell_to_scells = Table(tcell_to_scells_data,tcell_to_scells_ptr)
  k = CombineContributionsMap(scell_to_val)
  tcell_to_val = lazy_map(k,tcell_to_scells)
end

struct CombineContributionsMap{A} <: Map
  scell_to_val::A
end

function return_cache(k::CombineContributionsMap{<:AbstractArray{<:Number}},scells)
  z = zero(testitem(k.scell_to_val))
  z, array_cache(k.scell_to_val)
end

function evaluate!(cache,k::CombineContributionsMap{<:AbstractArray{<:Number}},scells)
  z,c = cache
  val = zero(z)
  for scell in scells
    val += getindex!(c,k.scell_to_val,scell)
  end
  val
end

function return_cache(k::CombineContributionsMap,scells)
  array = testitem(k.scell_to_val)
  c12 = _cache_compress(array)
  c = array_cache(k.scell_to_val)
  c12, c
end

function evaluate!(cache,k::CombineContributionsMap,scells)
  c12,c = cache
  c1,c2 = c12
  if length(scells) == 0
    _setempty_compress!(c1)
    val = _uncached_compress!(c1,c2)
  else
    val1 = getindex!(c,k.scell_to_val,scells[1])
    _setsize_compress!(c1,val1)
    val = _uncached_compress!(c1,c2)
    _copyto_compress!(val,val1)
    for i in 2:length(scells)
      vali = getindex!(c,k.scell_to_val,scells[i])
      _addto_compress!(val,vali)
    end
    val
  end
  val
end

function _cache_compress(array::AbstractArray)
  c1 = CachedArray(copy(array))
  c1, nothing
end

function _cache_compress(array::ArrayBlock)
  c1 = CachedArray(deepcopy(array))
  c2 = return_cache(Fields.unwrap_cached_array,c1)
  c1,c2
end

function _cache_compress(matvec::Tuple)
  mat, vec = matvec
  cmat,ca = _cache_compress(mat)
  cvec,cb = _cache_compress(vec)
  (cmat, cvec), (ca,cb)
end

function _setempty_compress!(a::CachedArray)
  setsize!(a,0 .* size(a))
end

function _setempty_compress!(a::ArrayBlock)
  for k in eachindex(a.array)
    if a.touched[k]
      _setempty_compress!(a.array[k])
    end
  end
end

function _setempty_compress!(ab::Tuple)
  a,b = ab
  _setempty_compress!(a)
  _setempty_compress!(b)
end

function _uncached_compress!(a::CachedArray,c2)
  a.array
end

function _uncached_compress!(c1::ArrayBlock,c2)
  evaluate!(c2,Fields.unwrap_cached_array,c1)
end

function _uncached_compress!(c1::Tuple,c2)
  a1,b1 = c1
  a2,b2 = c2
  a = _uncached_compress!(a1,a2)
  b = _uncached_compress!(b1,b2)
  (a,b)
end

function _setsize_compress!(a::CachedArray,b::AbstractArray)
  setsize!(a,size(b))
end

function _setsize_compress!(a::ArrayBlock,b::ArrayBlock)
  @check size(a) == size(b)
  for k in eachindex(a.array)
    @check a.touched[k] == b.touched[k]
    if a.touched[k]
    _setsize_compress!(a.array[k],b.array[k])
    end
  end
end

function _setsize_compress!(a::Tuple,b::Tuple)
  a1,a2 = a
  b1,b2 = b
  _setsize_compress!(a1,b1)
  _setsize_compress!(a2,b2)
end

function _copyto_compress!(a::AbstractArray,b::AbstractArray)
  @check size(a) == size(b)
  for k in eachindex(a)
    a[k] = b[k]
  end
  a
end

function _copyto_compress!(a::ArrayBlock,b::ArrayBlock)
  @check size(a) == size(b)
  for k in eachindex(a.array)
    @check a.touched[k] == b.touched[k]
    if a.touched[k]
      _copyto_compress!(a.array[k],b.array[k])
    end
  end
  a
end

function _copyto_compress!(a::Tuple,b::Tuple)
  a1,a2 = a
  b1,b2 = b
  _copyto_compress!(a1,b1)
  _copyto_compress!(a2,b2)
end

function _addto_compress!(a::AbstractArray,b::AbstractArray)
  @check size(a) == size(b)
  for k in eachindex(a)
    a[k] = a[k] + b[k]
  end
  a
end

function _addto_compress!(a::ArrayBlock,b::ArrayBlock)
  @check size(a) == size(b)
  for k in eachindex(a.array)
    @check a.touched[k] == b.touched[k]
    if a.touched[k]
      _addto_compress!(a.array[k],b.array[k])
    end
  end
  a
end

function _addto_compress!(a::Tuple,b::Tuple)
  a1,a2 = a
  b1,b2 = b
  _addto_compress!(a1,b1)
  _addto_compress!(a2,b2)
end


text:
@doc raw"""
    AffTransformation3

Represents three-dimensional affine transformations.

The general form of an affine transformation is based on a homogeneous
representation of points. Thereby all transformations can be realized by matrix
multiplications.

Multiplying the transformation matrix by a scalar does not change the
represented transformation. Therefore, any transformation represented by a
matrix with rational entries can be represented by a transformation matrix with
integer entries as well. (Multiply the matrix with the common denominator of
the rational entries.) Hence, it is sufficient to use the number type
[`RT`](@ref) to represent the entries of the transformation matrix.

CGAL offers several specialized affine transformations. Different constructors
are provided to create them. They are paremeterized with a symbolic name to
denote the transformation type, followed by additional parameters. The symbolic
name tags solve ambiguities in the function overloading and they make the code
more readable, i.e., what type of transformation is created.

In three-dimensional space we have a ``4 × 4`` matrix ``(m_{ij})_{i,j=0…3}``.
Entries ``m_{30}``, ``m_{31}``, and ``m_{30}`` are always zero and therefore do
not appear in the constructs.

See also: [`AffTransformation2`](@ref), [`IdentityTransformation`](@ref),
[`Rotation`](@ref), [`Scaling`](@ref), [`Translation`](@ref),
[`Reflection`](@ref), [`rational_rotation_approximation`](@ref)
"""
AffTransformation3

"""
    AffTransformation3(::IdentityTransformation)

Introduces an identity transformation.
"""
AffTransformation3(::IdentityTransformation)

"""
    AffTransformation3(::Translation, v::Vector3)

Introduces a translation by vector `v`.
"""
AffTransformation3(::Translation, v::Vector3)

"""
    AffTransformation3(::Scaling, s::Real, hw::Real = 1)

Introduces a scaling by a scale factor ``s/hw``.
"""
AffTransformation3(S::Scaling, s::Real, hw::Real = 1)

AffTransformation3(S, s::Real, hw::Real = 1) =
    AffTransformation3(S, convert(FT, s), convert(FT, hw))

@doc raw"""
    AffTransformation3(m₀₀::Real, m₀₁::Real, m₀₂::Real, m₀₃::Real,
                       m₁₀::Real, m₁₁::Real, m₁₂::Real, m₁₃::Real,
                       m₂₀::Real, m₂₁::Real, m₂₂::Real, m₂₃::Real,
                                                                   hw::Real = 1)

Introduces a general affine transformation in the matrix form ``
\small\left(\begin{array}{cccc}
    m_{00} & m_{01} & m_{02} & m_{03} \\ m_{10} & m_{11} & m_{12} & m_{13} \\ m_{20} & m_{21} & m_{22} & m_{23} \\ 0 & 0 & 0 & hw
\end{array}\right)
``.

The part ``
\frac{1}{hw}\small\left(\begin{array}{ccc}
    m_{00} & m_{01} & m_{02} \\ m_{10} & m_{11} & m_{12} \\ m_{20} & m_{21} & m_{22}
\end{array}\right)
``
defines the scaling and rotation part of the transformation, while the vector ``
\small\left(\begin{array}{c}
    m_{03} \\ m_{23} \\ m_{23}
\end{array}\right)
``
contains the translational part.
"""
AffTransformation3(m₀₀::Real, m₀₁::Real, m₀₂::Real, m₀₃::Real,
                   m₁₀::Real, m₁₁::Real, m₁₂::Real, m₁₃::Real,
                   m₂₀::Real, m₂₁::Real, m₂₂::Real, m₂₃::Real,
                                                               hw::Real = 1) =
    AffTransformation3(convert.(FT, (m₀₀, m₀₁, m₀₂, m₀₃,
                                     m₁₀, m₁₁, m₁₂, m₁₃,
                                     m₂₀, m₂₁, m₂₂, m₂₃,
                                                         hw))...)

@doc raw"""
    AffTransformation3(m₀₀::Real, m₀₁::Real, m₀₂::Real,
                       m₁₀::Real, m₁₁::Real, m₁₂::Real,
                       m₂₀::Real, m₂₁::Real, m₂₂::Real,
                                                        hw::Real = 1)

Introduces a general linear transformation ``
\small\left(\begin{array}{cccc}
    m_{00} & m_{01} & m_{02} & 0 \\ m_{10} & m_{11} & m_{12} & 0 \\ m_{20} & m_{21} & m_{22} & 0 \\ 0 & 0 & 0 & hw
\end{array}\right)
``,
i.e. an affine transformation without translational part.
"""
AffTransformation3(m₀₀::Real, m₀₁::Real, m₀₂::Real,
                   m₁₀::Real, m₁₁::Real, m₁₂::Real,
                   m₂₀::Real, m₂₁::Real, m₂₂::Real,
                                                    hw::Real = 1) =
    AffTransformation3(convert.(FT, (m₀₀, m₀₁, m₁₀, m₁₁, hw))...)

"""
    transform(t::AffTransformation3, x::Union{Point3,Vector3,Direction3,Line3})
    (::AffTransformation3)(x::Union{Point3,Vector3,Direction3,Line3})

The main thing to do with transformations is to apply them on geometric objects.

The transformations classes provide a [`transform`](@ref) function has methods
for points, vectors, directions, and planes. The same functionality is available
through the invocation of [`AffTransformation3`](@ref) instances.
"""
transform(t::AffTransformation3, p::Union{Point3,Vector3,Direction3,Plane3}),
(::AffTransformation3)(x::Union{Point3,Vector3,Direction3,Line3})

"""
    *(t::AffTransformation3, s::AffTransformation3)

Composes two affine transformations.
"""
*(t::AffTransformation3, s::AffTransformation3)

"""
    inverse(t::AffTransformation3)

Gives the inverse transformation.
"""
inverse(t::AffTransformation3)

"""
    ==(t::AffTransformation3, s::AffTransformation3)

Compares two affine transformations.
"""
==(t::AffTransformation3, s::AffTransformation3)

"""
    is_even(t::AffTransformation3)

Returns `true`, if the transformation is not reflecting, i.e., the determinant
of the involved linear transformation is non-negative.
"""
is_even(t::AffTransformation3)

"""
    is_odd(t::AffTransformation3)

Returns `true`, if the transformation is reflecting.
"""
is_odd(t::AffTransformation3)

"""
    is_scaling(t::AffTransformation3)

Returns `true`, if the object was constructed using the tag
[`CGAL.Scaling`](@ref), or is the result of the composition of only such scaling
transformation objects.
"""
is_scaling(t::AffTransformation3)

"""
    is_translation(t::AffTransformation3)

Returns `true`, if the object was constructed using the tag
[`CGAL.Translation`](@ref), or is the result of the composition of only such
translation transformation objects.
"""
is_translation(t::AffTransformation3)

"""
    m(t::AffTransformation3, i::Integer, j::Integer)

Returns entry ``m_{ij}`` in a matrix representation in which ``m_{22}`` is 1.
"""
m(t::AffTransformation3, i::Integer, j::Integer)

"""
    hm(t::AffTransformation3, i::Integer, j::Integer)

Returns entry ``m_{ij}`` in some fixed matrix representation.
"""
hm(t::AffTransformation3, i::Integer, j::Integer)


text:
function update_dmm_conf!(dmm_conf, dmm, refresh_cnt, base)
    # measurements
    refresh_cnt==base*1  && (dmm_conf.primary = string(get_primary_measurement(dmm, dmm_conf.crt_chan)))
    refresh_cnt==base*5  && begin
    if dmm_conf.crt_func == "RIPPLE"
        dmm_conf.secondary = string(get_secondary_measurement(dmm, dmm_conf.crt_chan))
    else
        dmm_conf.secondary = "----.---"
    end
    end
    # get current function
    refresh_cnt==base*10 && (dmm_conf.crt_func = get_sense_func(dmm, dmm_conf.crt_chan))
    # get range
    #refresh_cnt==base*15 && (dmm_conf.range = get_sense_range_auto(dmm, dmm_conf.crt_chan))
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
__precompile__(true)
module UnitfulMR

import Unitful
using Unitful: @unit, @u_str

# MR units
@unit Gauss "Gauss" Gauss 1e-4u"T" false # Gs occupied by Giga Sec in Unitful

# http://painterqubits.github.io/Unitful.jl/stable/extending/#Precompilation
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulMR)
end

end # module


text:
# some convenient helper function for Git CLI

function remote_repo_url(repo::String=pwd(); remote::String="origin")
    return cd(repo) do
        String(readchomp(`$(git()) config --get remote.$remote.url`))
    end
end

function remote_exists(repo::String=pwd())
    url = remote_repo_url(repo)
    p = cd(repo) do
        redirect_stdio(stderr=devnull, stdout=devnull) do
            run(ignorestatus(`$(git()) ls-remote --exit-code $url`))
        end
    end
    return p.exitcode == 0
end

function default_branch(repo::String=pwd(); remote::String="origin")
    remote_exists(repo) || error("cannot determine default branch, remote does not exist")

    url = remote_repo_url(repo)
    return cd(repo) do
        s = redirect_stdio(stderr=devnull) do
            readchomp(ignorestatus(`$(git()) ls-remote --symref $url HEAD`))
        end
        m = match(r"'s|^ref: refs/heads/(\S+)\s+HEAD|\1|p'", s)
        m === nothing && error("invalid remote response: $s")
        return String(m[1])
    end
end

function current_branch(repo::String=pwd())
    return cd(repo) do
        String(readchomp(`$(git()) rev-parse --abbrev-ref HEAD`))
    end
end

function checkout(repo::String=pwd(); branch::String, quiet::Bool=false)
    cd(repo) do
        if quiet
            run(`$(git()) checkout --quiet $branch`)
        else
            run(`$(git()) checkout $branch`)
        end
    end
end

function checkout(f, repo::String=pwd(); branch::String, quiet::Bool=false)
    current_br = current_branch(repo)
    checkout(repo; branch, quiet)
    ret = f()
    checkout(repo; branch=current_br, quiet)
    return ret
end

function git_clone(url::String, to::String)
    run(`$(git()) clone $url $to`)
end

function git_set_upstream(repo::String=pwd(); url::String)
    cd(repo) do
        run(`$(git()) remote add upstream $url`)
        run(`$(git()) fetch upstream`)
        main = default_branch(repo; remote="upstream")
        run(`$(git()) branch --set-upstream-to=upstream/$main`)
    end
end

function get_head_sha256(repo::String=pwd(); branch=current_branch(repo))
    return cd(repo) do
        String(readchomp(`$(git()) show -s --format="%H" $branch`))
    end
end

function list_branches(repo::String=pwd(); remote::Bool=false)
    return cd(repo) do
        results = readchomp(`$(git()) for-each-ref --format='%(refname:short)' refs/heads`)
        String.(split(results, '\n'))
    end
end

function git_get_toplevel_path(path::String=pwd())
    return cd(path) do
        String(readchomp(`git rev-parse --show-toplevel`))
    end
end

function isdiff(repo::String=pwd(); cached::Bool=false)
    return cd(repo) do
        if cached
            p = run(ignorestatus(`git diff --cached --quiet --exit-code`))
        else
            p = run(ignorestatus(`git diff --quiet --exit-code`))
        end
        return p.exitcode == 1
    end
end

isdirty(repo::String=pwd(); cached::Bool=false) = isdiff(repo; cached)

function github_repo(repo::String, remote="origin")
    url = remote_repo_url(repo; remote)
    github_https = "https://github.com/"
    github_ssh = "git@github.com:"
    if startswith(url, github_https)
        if endswith(url, ".git")
            return url[length(github_https)+1:end-4]
        else
            return url[length(github_https)+1:end]
        end
    elseif startswith(url, github_ssh)
        return url[length(github_ssh)+1:end-4]
    else
        return
    end
end

function fetch_repo_from_url(url; options...)
    HTTPS_GITHUB = "https://github.com/"
    GIT_GITHUB = "git@github.com:"
    if startswith(url, HTTPS_GITHUB) && endswith(url, ".git")
        repo = GitHub.repo(url[length(HTTPS_GITHUB)+1:end-4]; options...)
    elseif startswith(url, GIT_GITHUB) && endswith(url, ".git")
        repo = GitHub.repo(url[length(GIT_GITHUB)+1:end-4]; options...)
    else
        return "not a GitHub repo, please visit $(url) for details about this package"
    end

    return repo
end


text:
"""
functions that controll window - so basically treshords for mask display
"""
module WindowControll
using ModernGL, ..DisplayWords, ..StructsManag, Setfield, ..PrepareWindow,   ..DataStructs , Rocket, GLFW,Dictionaries,  ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils,  ..Uniforms, Match, Parameters,DataTypesBasic
export  setTextureWindow,getNewTresholdValue

"""
KEY_F1 - will display wide window for bone Int32(1000),Int32(-1000)
KEY_F2 - will display window for soft tissues Int32(400),Int32(-200)
KEY_F3 - will display wide window for lung viewing  Int32(0),Int32(-1000)
KEY_F4,  KEY_F5 -
    sets minimum (F4) and maximum (KEY_F5) value for display (with combination of + and minus signs - to increase or decrease given treshold) -
        in case of continuus colors it will clamp values - so all above max will be equaled to max ; and min if smallert than min
        in case of main CT mask - it will controll min shown white and max shown black
        in case of maks with single color associated we will step data so if data is outside the rande it will return 0 - so will not affect display
KEY_F6 - controlls contribution  of given mask to the overall image - maximum value is 1 minimum 0 if we have 3 masks and all control contribution is set to 1 and all are visible their corresponding influence to pixel color is 33%
      if plus is pressed it will increse contribution by 0.1
      if minus is pressed it will decrease contribution by 0.1
"""
function processKeysInfo(wind::Identity{WindowControlStruct}
    ,actor::SyncActor{Any, ActorWithOpenGlObjects}
    ,keyInfo::KeyboardStruct
    ,toBeSavedForBack::Bool = true) where T
    #we have some predefined windows
    joined = join(keyInfo.lastKeysPressed)

    old = actor.actor.mainForDisplayObjects.windowControlStruct
    windowStruct =  primaryModificationsOfWindContr(wind.value,keyInfo)

    dispatchToFunctions(windowStruct,actor,keyInfo)

    #to display change

   basicRender(actor.actor.mainForDisplayObjects.window)

       # for undoing action
    if(toBeSavedForBack)
         addToforUndoVector(actor, ()-> processKeysInfo( Option(old),actor, keyInfo,false ))
    end

end#processKeysInfo

"""
On the basis of the input WindowControlStruct and keyInfo it makes necessary primary modifications to WindowControlStruct
"""
function primaryModificationsOfWindContr(windowStruct::WindowControlStruct,keyInfo::KeyboardStruct )::WindowControlStruct
    return @match windowStruct.letterCode begin
    "F1" => WindowControlStruct(letterCode="F1",min_shown_white=Int32(1000), max_shown_black=Int32(-1000) )
    "F2" => WindowControlStruct(letterCode="F2",min_shown_white=Int32(400), max_shown_black=Int32(-200) )
    "F3" => WindowControlStruct(letterCode="F3",min_shown_white=Int32(0),max_shown_black= Int32(-1000) )
    "F4" => WindowControlStruct(letterCode="F4", toIncrease= keyInfo.isPlusPressed, toDecrease= keyInfo.isMinusPressed,lower=true )
    "F5" => WindowControlStruct(letterCode="F5",  toIncrease= keyInfo.isPlusPressed,toDecrease=keyInfo.isMinusPressed,upper=true)
    "F6" => WindowControlStruct(letterCode="F6",  toIncrease= keyInfo.isPlusPressed,toDecrease=keyInfo.isMinusPressed,maskContributionToChange=true)
    _ => windowStruct
        end

end   #primaryModificationsOfWindContr

"""
Based on window struct and key info it will controll  which function should be invoked

"""
function dispatchToFunctions(windowStruct::WindowControlStruct
    ,actor::SyncActor{Any, ActorWithOpenGlObjects}
    ,keyInfo::KeyboardStruct )

   mainWindows=  @match (windowStruct.letterCode) begin
        "F1" =>  (setmainWindow(actor,windowStruct), "Fsth")
        "F2" =>  (setmainWindow(actor,windowStruct), "Fsth")
        "F3" =>  (setmainWindow(actor,windowStruct), "Fsth")
        bad  => "nothing"
end#match

 textureList= actor.actor.textureToModifyVec
currentWindowInActor = actor.actor.mainForDisplayObjects.windowControlStruct
if(mainWindows== "nothing" && !isempty(textureList))
    textur = textureList[1]
     matched = @match (windowStruct.letterCode, windowStruct.toIncrease ,  windowStruct.toDecrease ) begin
        ("F4" ,true,false) =>    lowTreshUp(windowStruct,actor,textur,currentWindowInActor)
        ("F4" ,false,true) =>    lowTreshDown(windowStruct,actor,textur,currentWindowInActor)
        ("F5" ,false,true) =>    highTreshDown(windowStruct,actor,textur,currentWindowInActor)
        ("F5" ,true,false) =>    highTreshUp(windowStruct,actor,textur,currentWindowInActor)
        ("F6" ,true,false) =>    maskContrUp(windowStruct,actor,textur,currentWindowInActor)
        ("F6" ,false,true) =>    maskContrDown(windowStruct,actor,textur,currentWindowInActor)
        bad                 => "nothing"
    end#match
else
        #updating current windowing object and getting reference to old
        actor.actor.mainForDisplayObjects= setproperties(actor.actor.mainForDisplayObjects,(windowControlStruct =windowStruct))

end #if

end#dispatchToFunctions

"""
sets lower treshold and Increase it
"""
function lowTreshUp(windowStruct::WindowControlStruct,actor::SyncActor{Any, ActorWithOpenGlObjects},textur::TextureSpec,currentWindowInActor::WindowControlStruct)
    if(textur.isMainImage)
        setmainWindow(actor,setproperties(currentWindowInActor, (max_shown_black= currentWindowInActor.max_shown_black+15)))
    else
        textur.minAndMaxValue[1]+=getNewTresholdChangeValue(textur)
        setTextureWindow( textur, actor )
    end#if
end#lowTreshUp
"""
sets lower treshold and decrese it
"""
function lowTreshDown(windowStruct::WindowControlStruct,actor::SyncActor{Any, ActorWithOpenGlObjects},textur::TextureSpec,currentWindowInActor::WindowControlStruct)
    if(textur.isMainImage)
        setmainWindow(actor,setproperties(currentWindowInActor, (max_shown_black= currentWindowInActor.max_shown_black-15)))
    else
        texturParam = parameter_type(textur)
        if(texturParam==UInt8 ||texturParam==UInt16 ||texturParam==UInt32 ||texturParam==UInt64 )
            textur.minAndMaxValue[1]= max(0,textur.minAndMaxValue[1]-getNewTresholdChangeValue(textur))
        else
            textur.minAndMaxValue[1]-=getNewTresholdChangeValue(textur)
        end#if
        setTextureWindow( textur, actor )
    end#if
end#lowTreshDown

"""
sets upper treshold and Increase it
"""
function highTreshUp(windowStruct::WindowControlStruct,actor::SyncActor{Any, ActorWithOpenGlObjects},textur::TextureSpec,currentWindowInActor::WindowControlStruct)
    if(textur.isMainImage)
        setmainWindow(actor,setproperties(currentWindowInActor, (min_shown_white= currentWindowInActor.min_shown_white+15)))
    else

            textur.minAndMaxValue[2]+=max(getNewTresholdChangeValue(textur),1)

        setTextureWindow( textur, actor )
    end#if
end#highTreshUp

"""
sets upper treshold and decrese it
"""
function highTreshDown(windowStruct::WindowControlStruct,actor::SyncActor{Any, ActorWithOpenGlObjects},textur::TextureSpec,currentWindowInActor::WindowControlStruct)
    if(textur.isMainImage)
        setmainWindow(actor,setproperties(currentWindowInActor, (min_shown_white= currentWindowInActor.min_shown_white-15)))
    else
        texturParam = parameter_type(textur)
        if(texturParam==UInt8 ||texturParam==UInt16 ||texturParam==UInt32 ||texturParam==UInt64 )
            textur.minAndMaxValue[2]= maximum([0,textur.minAndMaxValue[2]-getNewTresholdChangeValue(textur), textur.minAndMaxValue[1]] )
        else
            textur.minAndMaxValue[2]=max(textur.minAndMaxValue[2]-getNewTresholdChangeValue(textur) ,textur.minAndMaxValue[1] )
        end#if
        setTextureWindow( textur, actor )
    end#if
end#highTreshDown

"""
sets mask contribution and  decrese it
"""
function maskContrDown(windowStruct::WindowControlStruct,actor::SyncActor{Any, ActorWithOpenGlObjects},textur::TextureSpec,currentWindowInActor::WindowControlStruct)
    if(textur.isMainImage)
        changeMainTextureContribution(textur,Float32(-0.1), actor)
    else
        changeTextureContribution(textur,Float32(-0.1) )
    end#if
end#maskContrDown

"""
sets mask contribution and increase it
"""
function maskContrUp(windowStruct::WindowControlStruct,actor::SyncActor{Any, ActorWithOpenGlObjects},textur::TextureSpec,currentWindowInActor::WindowControlStruct)
    if(textur.isMainImage)
        changeMainTextureContribution(textur, Float32(0.1), actor)
    else
        changeTextureContribution(textur, Float32(0.1) )
    end#if
end#maskContrUp

"""
set main window - min shown white and max shown black on the basis of textur data and windowStruct
"""
function setmainWindow(actor::SyncActor{Any, ActorWithOpenGlObjects},windowStruct::WindowControlStruct)
        #updating current windowing object and getting reference to old
        setCTWindow(windowStruct.min_shown_white,windowStruct.max_shown_black, actor.actor.mainForDisplayObjects.mainImageUniforms)
        actor.actor.mainForDisplayObjects= setproperties(actor.actor.mainForDisplayObjects,(windowControlStruct =windowStruct))
end#setmainWindow

"""
sets minimum and maximum value for display -
    in case of continuus colors it will clamp values - so all above max will be equaled to max ; and min if smallert than min
    in case of main CT mask - it will controll min shown white and max shown black
    in case of maks with single color associated we will step data so if data is outside the rande it will return 0 - so will not affect display
"""
function setTextureWindow(textur::TextureSpec,actor::SyncActor{Any, ActorWithOpenGlObjects})
    # activeTextureList= actor.actor.textureToModifyVec
    # allTexturesList = actor.actor.mainForDisplayObjects.listOfTextSpecifications
    # notModifiedTextures= filter(it->it.name!=textur.name, allTexturesList)

    coontrolMinMaxUniformVals(textur)
end#    setTextureWindow

"""
helper function for setTextureWindow on the basis of given texture spec will give value proportional to the range
"""
function getNewTresholdChangeValue(textur::TextureSpec)::Int64
    return  ceil(abs(textur.minAndMaxValue[2]-textur.minAndMaxValue[1])/20 )
end #getNewTresholdValue

end#WindowControll

    # #in case we pressed F4 or F5 we want to manually change the window
    #     if( !isempty(textureList) && (windowStruct.toIncrease || windowStruct.toDecrease    ) )
    #         textur= textureList[1]
    #         if(textur.isMainImage)
    #             if(windStr.lower)
    #                 windowStruct= setproperties(windowStruct,  (max_shown_black= getNewTresholdValue(old.max_shown_black, windowStruct),min_shown_white=old.min_shown_white    ) )
    #             elseif(windStr.upper)
    #                 windowStruct= setproperties(windowStruct,  (min_shown_white= getNewTresholdValue(old.min_shown_white, windowStruct) , max_shown_black=old.max_shown_black  ) )
    #             end#ifs
    #             setCTWindow(windowStruct.min_shown_white,windowStruct.max_shown_black, actor.actor.mainForDisplayObjects.mainImageUniforms)
    #         else# if texture is not main image
    #             setTextureWindow(textur, windowStruct)
    #         end#if
    #     else
    #     # setting window and showig it
    #     setCTWindow(windowStruct.min_shown_white,windowStruct.max_shown_black, actor.actor.mainForDisplayObjects.mainImageUniforms)
    #     end#if


text:
# # Dry atmosphere GCM with Held-Suarez forcing
#
# The Held-Suarez setup (Held and Suarez, 1994) is a textbook example for a
# simplified atmospheric global circulation model configuration which has been
# used as a benchmark experiment for development of the dynamical cores (i.e.,
# GCMs without continents, moisture or parametrization schemes of the physics)
# for atmospheric models.  It is forced by a thermal relaxation to a reference
# state and damped by linear (Rayleigh) friction. This example demonstrates how
#
#   * to set up a ClimateMachine-Atmos GCM configuration;
#   * to select and save GCM diagnostics output.
#
# To begin, we load ClimateMachine and a few miscellaneous useful Julia packages.
using Distributions
using Random
using StaticArrays
using UnPack

# ClimateMachine specific modules needed to make this example work (e.g., we will need
# spectral filters, etc.).
using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.TemperatureProfiles
using ClimateMachine.SystemSolvers
using ClimateMachine.ODESolvers
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using ClimateMachine.BalanceLaws

import ClimateMachine.BalanceLaws: source, prognostic_vars

# [ClimateMachine parameters](https://github.com/CliMA/CLIMAParameters.jl) are
# needed to have access to Earth's physical parameters.
using CLIMAParameters
using CLIMAParameters.Planet: MSLP, R_d, day, grav, cp_d, cv_d, planet_radius

# We need to load the physical parameters for Earth to have an Earth-like simulation :).
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet();

"""
    HeldSuarezForcingTutorial <: TendencyDef{Source}

Defines a forcing that parametrises radiative and frictional effects using
Newtonian relaxation and Rayleigh friction, following Held and Suarez (1994)
"""
struct HeldSuarezForcingTutorial <: TendencyDef{Source} end

prognostic_vars(::HeldSuarezForcingTutorial) = (Momentum(), Energy())

function held_suarez_forcing_coefficients(bl, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed
    FT = eltype(state)

    ## Parameters
    T_ref = FT(255)
    param_set = parameter_set(bl)

    _R_d = FT(R_d(param_set))
    _day = FT(day(param_set))
    _grav = FT(grav(param_set))
    _cp_d = FT(cp_d(param_set))
    _p0 = FT(MSLP(param_set))

    ## Held-Suarez parameters
    k_a = FT(1 / (40 * _day))
    k_f = FT(1 / _day)
    k_s = FT(1 / (4 * _day))
    ΔT_y = FT(60)
    Δθ_z = FT(10)
    T_equator = FT(315)
    T_min = FT(200)
    σ_b = FT(7 / 10)

    ## Held-Suarez forcing
    φ = latitude(bl, aux)
    p = air_pressure(ts)

    ##TODO: replace _p0 with dynamic surface pressure in Δσ calculations to account
    ##for topography, but leave unchanged for calculations of σ involved in T_equil
    σ = p / _p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor
    return (k_v = k_v, k_T = k_T, T_equil = T_equil)
end

function source(::Energy, s::HeldSuarezForcingTutorial, atmos, args)
    @unpack state = args
    FT = eltype(state)
    @unpack ts = args.precomputed
    nt = held_suarez_forcing_coefficients(atmos, args)
    param_set = parameter_set(atmos)
    _cv_d = FT(cv_d(param_set))
    @unpack k_T, T_equil = nt
    T = air_temperature(ts)
    return -k_T * state.ρ * _cv_d * (T - T_equil)
end

function source(::Momentum, s::HeldSuarezForcingTutorial, atmos, args)
    @unpack state, aux = args
    nt = held_suarez_forcing_coefficients(atmos, args)
    return -nt.k_v * projection_tangential(atmos, aux, state.ρu)
end

# ## Set initial condition
# When using ClimateMachine, we need to define a function that sets the initial
# state of our model run. In our case, we use the reference state of the
# simulation (defined below) and add a little bit of noise. Note that the
# initial states includes a zero initial velocity field.
function init_heldsuarez!(problem, balance_law, state, aux, localgeo, time)
    FT = eltype(state)

    ## Set initial state to reference state with random perturbation
    rnd = FT(1 + rand(Uniform(-1e-3, 1e-3)))
    state.ρ = aux.ref_state.ρ
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.energy.ρe = rnd * aux.ref_state.ρe
end;

# ## Initialize ClimateMachine
# Before we do anything further, we need to initialize ClimateMachine. Among
# other things, this will initialize the MPI us.
ClimateMachine.init();

# ## Setting the floating-type precision
# ClimateMachine allows us to run a model with different floating-type
# precisions, with lower precision we get our results faster, and with higher
# precision, we may get more accurate results, depending on the questions we
# are after.
const FT = Float64;

# ## Setup model configuration
# Now that we have defined our forcing and initialization functions, and have
# initialized ClimateMachine, we can set up the model.
#
# ## Set up a reference state
# We start by setting up a reference state. This is simply a vector field that
# we subtract from the solutions to the governing equations to both improve
# numerical stability of the implicit time stepper and enable faster model
# spin-up. The reference state assumes hydrostatic balance and ideal gas law,
# with a pressure $p_r(z)$ and density $\rho_r(z)$ that only depend on altitude
# $z$ and are in hydrostatic balance with each other.
#
# In this example, the reference temperature field smoothly transitions from a
# linearly decaying profile near the surface to a constant temperature profile
# at the top of the domain.
temp_profile_ref = DecayingTemperatureProfile{FT}(param_set)
ref_state = HydrostaticState(temp_profile_ref);

# ## Set up a Rayleigh sponge layer
# To avoid wave reflection at the top of the domain, the model applies a sponge
# layer that linearly damps the momentum equations.
domain_height = FT(30e3)               # height of the computational domain (m)
z_sponge = FT(12e3)                    # height at which sponge begins (m)
α_relax = FT(1 / 60 / 15)              # sponge relaxation rate (1/s)
exponent = FT(2)                       # sponge exponent for squared-sinusoid profile
u_relax = SVector(FT(0), FT(0), FT(0)) # relaxation velocity (m/s)
sponge = RayleighSponge{FT}(domain_height, z_sponge, α_relax, u_relax, exponent);

# ## Set up turbulence models
# In order to produce a stable simulation, we need to dissipate energy and
# enstrophy at the smallest scales of the developed flow field. To achieve this
# we set up diffusive forcing functions.
c_smag = FT(0.21);   # Smagorinsky constant
τ_hyper = FT(4 * 3600); # hyperdiffusion time scale
turbulence_model = SmagorinskyLilly(c_smag);
hyperdiffusion_model = DryBiharmonic(FT(4 * 3600));

# ## Instantiate the model
# The Held Suarez setup was designed to produce an equilibrated state that is
# comparable to the zonal mean of the Earth’s atmosphere.
physics = AtmosPhysics{FT}(
    param_set;
    ref_state = ref_state,
    turbulence = turbulence_model,
    hyperdiffusion = hyperdiffusion_model,
    moisture = DryModel(),
);

model = AtmosModel{FT}(
    AtmosGCMConfigType,
    physics;
    init_state_prognostic = init_heldsuarez!,
    source = (Gravity(), Coriolis(), HeldSuarezForcingTutorial(), sponge),
);

# This concludes the setup of the physical model!

# ## Set up the driver
# We just need to set up a few parameters that define the resolution of the
# discontinuous Galerkin method and for how long we want to run our model
# setup.
poly_order = 5;                        ## discontinuous Galerkin polynomial order
n_horz = 2;                            ## horizontal element number
n_vert = 2;                            ## vertical element number
resolution = (n_horz, n_vert)
n_days = 0.1;                          ## experiment day number
timestart = FT(0);                     ## start time (s)
timeend = FT(n_days * day(param_set)); ## end time (s);

# The next lines set up the spatial grid.
driver_config = ClimateMachine.AtmosGCMConfiguration(
    "HeldSuarez",
    poly_order,
    resolution,
    domain_height,
    param_set,
    init_heldsuarez!;
    model = model,
);

# The next lines set up the time stepper. Since the resolution
# in the vertical is much finer than in the horizontal,
# the 'stiff' parts of the PDE will be in the vertical.
# Setting `splitting_type = HEVISplitting()` will treat
# vertical acoustic waves implicitly, while all other dynamics
# are treated explicitly.
ode_solver_type = ClimateMachine.IMEXSolverType(
    splitting_type = HEVISplitting(),
    implicit_model = AtmosAcousticGravityLinearModel,
    implicit_solver = ManyColumnLU,
    solver_method = ARK2GiraldoKellyConstantinescu,
);

solver_config = ClimateMachine.SolverConfiguration(
    timestart,
    timeend,
    driver_config,
    Courant_number = FT(0.1),
    ode_solver_type = ode_solver_type,
    init_on_cpu = true,
    CFL_direction = HorizontalDirection(),
    diffdir = HorizontalDirection(),
);

# ## Set up spectral exponential filter
# After every completed time step we apply a spectral filter to remove
# remaining small-scale noise introduced by the numerical procedures. This
# assures that our run remains stable.
filterorder = 10;
filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder);
cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
    Filters.apply!(
        solver_config.Q,
        AtmosFilterPerturbations(model),
        solver_config.dg.grid,
        filter,
        state_auxiliary = solver_config.dg.state_auxiliary,
    )
    nothing
end;

# ## Setup diagnostic output
#
# Choose frequency and resolution of output, and a diagnostics group (dgngrp)
# which defines output variables. This needs to be defined in
# [`Diagnostics`](@ref ClimateMachine.Diagnostics).
interval = "1000steps";
_planet_radius = FT(planet_radius(param_set));
info = driver_config.config_info;
boundaries = [
    FT(-90.0) FT(-180.0) _planet_radius
    FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
];
resolution = (FT(10), FT(10), FT(1000)); # in (deg, deg, m)
interpol = ClimateMachine.InterpolationConfiguration(
    driver_config,
    boundaries,
    resolution,
);

dgngrps = [
    setup_dump_state_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    ),
    setup_dump_aux_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    ),
];
dgn_config = ClimateMachine.DiagnosticsConfiguration(dgngrps);

# ## Run the model
# Finally, we can run the model using the physical setup and solvers from
# above. We use the spectral filter in our callbacks after every time step, and
# collect the diagnostics output.
result = ClimateMachine.invoke!(
    solver_config;
    diagnostics_config = dgn_config,
    user_callbacks = (cbfilter,),
    check_euclidean_distance = true,
);

# ## References
#
# - [Held1994](@cite)


text:
using SeisIO, JLD2, SeisNoise, Statistics, PlotlyJS
include("utils.jl")
include("io.jl")
include("stacking.jl")

"""
    compute_reference_xcorr(tstamp::String, basefiname::String; phase_smoothing::Float64=0., stack::String="linear", reference::String="", thresh::Float64=-1.0)

Stack all cross-correlation functions for all given station pairs to generate a reference cross-correlation.

# Arguments
- `tstamp::String`    :
- `basefiname::String,`    : Input base file name e.g. "./inputData/BPnetwork"
- `phase_smoothing::Float64`    : Level of phase_smoothing (0 for linear stacking)
- `stack::String`     : "selective" for selective stacking and "linear" for linear stacking
- `reference::String`    : Path to the reference used in selective stacking
- `thresh::Float64`     : Threshold used for selective stacking

# Output
- `foname.jld2`    : contains arrays of reference cross-correlations for each station pair
"""
function compute_reference_xcorr(tstamp::String, basefiname::String; phase_smoothing::Float64=0., stack::String="linear", reference::String="", thresh::Float64=-1.0)
    # hold reference xcorrs in memory and write all at once
    ref_dict = Dict()

    # read unstacked xcorrs for each time stamp
    f_cur = jldopen(basefiname*".$tstamp.jld2")
    grp = f_cur[tstamp] # xcorrs
    println("$tstamp")

    # iterate over station pairs
    for pair in sort(keys(grp))
        # TODO: use only unique station pairings when creating references. Currently no guarantee of uniqueness (reverse can exist)
        # load xcorr
        xcorr = try grp[pair] catch; continue end
        #remove_nan!(xcorr)

        # stack xcorrs over length of CorrData object using either "selective" stacking or "linear" stacking
        if stack=="selective"
            f_exref = jldopen(reference)
            ref = f_exref[pair]
            close(f_exref)
            xcorr, rmList = selective_stacking(xcorr, ref, threshold=thresh)
        else
            stack!(xcorr, allstack=true, phase_smoothing=phase_smoothing)
        end

        # stack xcorrs if they have a key, assign key if not
        if haskey(ref_dict, pair)
            ref_dict[pair].corr .+= xcorr.corr
        else
            ref_dict[pair] = xcorr
        end
    end
    close(f_cur) # current xcorr file

    return ref_dict
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
@time @testset "BufferedConfigRecorder" begin
    @testset "Runtime File IO" begin
        isfile("testfile.confs") && rm("testfile.confs")
        r = MonteCarlo.BufferedConfigRecorder{Matrix{Float64}}("testfile.confs", 2, 20)

        # Initialization
        @test r.filename == MonteCarlo.FilePath(false, "testfile.confs", "testfile.confs")
        @test length(r.buffer) == 20
        @test eltype(r.buffer) == Matrix{Float64}
        @test r.rate == 2
        @test r.idx == 1
        @test r.chunk == 1
        @test r.total_length == 0
        @test r.save_idx == -1
        @test !isdefined(r.buffer, 1)

        # Basic utility functions
        @test length(r) == r.total_length
        @test lastindex(r) == r.total_length
        @test isempty(r)

        # push to empty memory buffer
        c = rand(4, 4)
        MonteCarlo._push!(r, c)

        @test r.idx == 2
        @test r.chunk == 1
        @test r.total_length == 1
        @test r.save_idx == -1
        @test r.buffer[1] == c
        @test !isdefined(r.buffer, 2)

        @test length(r) == r.total_length
        @test lastindex(r) == r.total_length
        @test !isempty(r)
        @test r[1] == c

        # fill memory buffer - no save yet
        for _ in 2:20
            MonteCarlo._push!(r, c)
        end

        @test r.idx == 21
        @test r.chunk == 1
        @test r.total_length == 20
        @test r.save_idx == -1
        @test r.buffer[20] == c
        @test r[20] == c

        # overflow push - should save and reset buffer
        c2 = rand(4, 4)
        MonteCarlo._push!(r, c2)

        @test r.idx == 2
        @test r.chunk == 2
        @test r.total_length == 21
        @test r.save_idx == 20
        @test r.buffer[1] == c2
        @test r[21] == c2

        # query last chunk - should load and mark idx for reload
        @test r[20] == c
        @test r[1] == c
        @test r.idx == -1
        @test r.chunk == 1
        @test r.total_length == 21
        @test r.save_idx == 21

        # push to wrong chunk - should load correct chunk
        c3 = rand(4, 4)
        MonteCarlo._push!(r, c3)

        @test r.idx == 3
        @test r.chunk == 2
        @test r.total_length == 22
        @test r.save_idx == 21
        @test r.buffer[1] == c2 # this would fail if chunk not loaded
        @test r.buffer[2] == c3
        @test r.buffer[3] == c

        # test merge!
        c4 = rand(4, 4)
        cs = ConfigRecorder{Matrix{Float64}}([c4 for _ in 1:17], 10)
        merge!(r, cs)

        @test r.idx == 20
        @test r.chunk == 2
        @test r.total_length == 39
        @test r.save_idx == 21
        @test r.buffer[1] == c2   # old
        @test r.buffer[2] == c3   # old
        @test r.buffer[3] == c4   # new
        @test r.buffer[19] == c4  # new
        @test r.buffer[20] == c   # to be overwritten
        rm("testfile.confs")
    end

    dir1 = randstring(16)
    dir2 = randstring(16)
    mkdir(dir1)
    mkdir(dir2)

    @testset "Parent related file IO (moving, renaming, replacing)" begin
        m = HubbardModel(2, 2, U = 1.0)
        recorder = BufferedConfigRecorder(DensityHirschField, RelativePath("testfile.confs"))
        link_id = recorder.link_id
        mc = DQMC(m, beta=1.0, recorder = recorder)
        mc.field.conf .= rand(mc.field)
        push!(recorder, mc.field, 0)

        # BCR should not save on construction
        @test recorder.filename.absolute_path == joinpath(pwd(), "testfile.confs")
        @test !isfile("testfile.confs")

        # parent save should trigger save
        MonteCarlo.save("$dir1/testfile.jld2", mc)
        @test isfile("$dir1/testfile.confs")
        @test isfile("$dir1/testfile.jld2")
        @test !isfile("testfile.confs")

        # both files should have same link_id
        MonteCarlo.JLD2.jldopen("$dir1/testfile.confs", "r") do file
            @test file["link_id"] == link_id
        end
        MonteCarlo.JLD2.jldopen("$dir1/testfile.jld2", "r") do file
            @test file["MC"]["configs"]["link_id"] == link_id
        end

        # check if data loaded is correct
        x = MonteCarlo.load("$dir1/testfile.jld2")
        @test x.recorder.filename == recorder.filename
        @test x.recorder.link_id == link_id
        @test recorder[1] == x.recorder[1]

        function move_load_check(recorder, filename_should_match, source_exists = false)
            mv("$dir1/testfile.jld2", "$dir2/testfile.jld2")
            x = MonteCarlo.load("$dir2/testfile.jld2")
            @test filename_should_match == (x.recorder.filename.relative_path == recorder.filename.relative_path)
            @test x.recorder.filename.absolute_path == joinpath(pwd(), dir2, x.recorder.filename.relative_path)
            @test x.recorder.link_id == link_id
            @test recorder[1] == x.recorder[1]
            filename = x.recorder.filename.relative_path
            @test isfile("$dir2/$filename")
            @test source_exists == isfile("$dir1/testfile.confs")
            @test length(readdir(dir2)) == 3 - filename_should_match
            mv("$dir2/testfile.jld2", "$dir1/testfile.jld2", force = true)
            mv("$dir2/$filename", "$dir1/testfile.confs", force = true)
            nothing
        end

        # Move parent only - recorder should move on load
        move_load_check(recorder, true)

        # Move both - recorder should adjust filename
        mv("$dir1/testfile.confs", "$dir2/testfile.confs")
        move_load_check(recorder, true)

        # Move parent, dublicate BCS file
        # file should be recognized as the same and both should stay
        cp("$dir1/testfile.confs", "$dir2/testfile.confs")
        move_load_check(recorder, true, true)

        # Move parent, create collision
        # file should be renamed
        open(f -> write(f, "Test"), "$dir2/testfile.confs", "w")
        try
            move_load_check(recorder, false)
        catch e
            @warn "Windows please" exception = e
        end

        # save uses the same function to move/replace/rename so it's not really
        # necessary to test this specifically
    end

    try
        rm(dir1, recursive = true, force = true)
        rm(dir2, recursive = true, force = true)
    catch e
        @warn "Windows please" exception = e
    end
end

function test_mc(mc, x)
    # Check if loaded/replayed mc matches original
    for f in fieldnames(typeof(mc.p))
        @test getfield(mc.p, f) == getfield(x.p, f)
    end
    @test mc.conf == x.conf
    @test mc.model.L == x.model.L
    @test mc.model.dims == x.model.dims
    for f in fieldnames(typeof(mc.model.l))
        @test getfield(mc.model.l, f) == getfield(x.model.l, f)
    end
    # @test mc.model.neighs == x.model.neighs
    @test mc.model.energy[] == x.model.energy[]
    for (k, v) in mc.thermalization_measurements
        for f in fieldnames(typeof(v))
            r = getfield(v, f) == getfield(x.thermalization_measurements[k], f)
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    for (k, v) in mc.measurements
        for f in fieldnames(typeof(v))
            r = getfield(v, f) == getfield(x.measurements[k], f)
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    nothing
end

@time @testset "MC" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(
        model, beta = 0.66, thermalization = 33, sweeps = 123,
        recorder = ConfigRecorder(MC, IsingModel, 1)
    )
    run!(mc, verbose=false)
    save("testfile.jld2", mc)
    x = load("testfile.jld2")
    rm("testfile.jld2")
    test_mc(mc, x)

    x.measurements = MonteCarlo.default_measurements(mc, model)
    x.last_sweep = 0
    replay!(x, verbose=false)
    test_mc(mc, x)

    # Test resume

    # Run for 1s with known RNG
    Random.seed!(123)
    model = IsingModel(dims=2, L=10)
    mc = MC(
        model, beta=1.0, thermalization = 0, sweeps=10_000_000,
        measure_rate=10_000, recorder=ConfigRecorder(MC, IsingModel, 10_000)
    )
    state = run!(
        mc, verbose = false,
        safe_before = now() + Second(1),
        grace_period = Millisecond(0),
        resumable_filename = "resumable_testfile.jld"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) > 1 "No measurements have been taken. Test with more time!"
    L = length(cs)

    # Test whether safe file gets overwritten correctly
    mc, state = resume!(
        "resumable_testfile.jld",
        verbose = false,
        safe_before = now() + Second(10),
        grace_period = Millisecond(0),
        overwrite = true,
        resumable_filename = "resumable_testfile.jld"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) - L > 1 "No new measurements have been taken. Test with more time!"
    @test isfile("resumable_testfile.jld")

    # Test whether data from resumed simulation is correct
    Random.seed!(123)
    model = IsingModel(dims=2, L=10)
    mc = MC(
        model, beta=1.0,
        thermalization=0, sweeps=10_000length(cs),
        measure_rate=10_000, recorder=ConfigRecorder(MC, IsingModel, 10_000)
    )
    state = run!(mc, verbose = false)
    @test mc.configs.configs == cs.configs
    @test mc.configs.rate == cs.rate
    rm("resumable_testfile.jld")
end

function test_dqmc(mc, x)
    for f in fieldnames(typeof(mc.parameters))
        @test getfield(mc.parameters, f) == getfield(x.parameters, f)
    end
    @test mc.field.conf == x.field.conf
    @test mc.model.mu == x.model.mu
    @test mc.model.t == x.model.t
    @test mc.model.U == x.model.U
    for f in fieldnames(typeof(mc.model.l))
        @test getfield(mc.model.l, f) == getfield(x.model.l, f)
    end
    @test MonteCarlo.nflavors(mc.field) == MonteCarlo.nflavors(x.field)
    @test mc.scheduler == x.scheduler
    for (k, v) in mc.thermalization_measurements
        for f in fieldnames(typeof(v))
            r = if getfield(v, f) isa LightObservable
                # TODO
                # implement == for LightObservable in MonteCarloObservable
                getfield(v, f).B == getfield(x.measurements[k], f).B
            else
                getfield(v, f) == getfield(x.measurements[k], f)
            end
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    for (k, v) in mc.measurements
        for f in fieldnames(typeof(v))
            v isa MonteCarlo.DQMCMeasurement && f == :temp && continue
            v isa MonteCarlo.DQMCMeasurement && f == :kernel && continue
            r = if getfield(v, f) isa LightObservable
                # TODO
                # implement == for LightObservable in MonteCarloObservable
                # TODO: implement ≈ for LightObservable, LogBinner, etc
                r = true
                a = getfield(v, f)
                b = getfield(x.measurements[k], f)
                for i in eachindex(getfield(v, f).B.compressors)
                    r = r && (a.B.compressors[i].value ≈ b.B.compressors[i].value)
                    r = r && (a.B.compressors[i].switch ≈ b.B.compressors[i].switch)
                end
                r = r && (a.B.x_sum ≈ b.B.x_sum)
                r = r && (a.B.x2_sum ≈ b.B.x2_sum)
                r = r && (a.B.count ≈ b.B.count)
            elseif getfield(v, f) isa LogBinner
                r = true
                a = getfield(v, f)
                b = getfield(x.measurements[k], f)
                for i in eachindex(a.compressors)
                    r = r && (a.compressors[i].value ≈ b.compressors[i].value)
                    r = r && (a.compressors[i].switch ≈ b.compressors[i].switch)
                end
                r = r && (a.x_sum ≈ b.x_sum)
                r = r && (a.x2_sum ≈ b.x2_sum)
                r = r && (a.count ≈ b.count)
            else
                getfield(v, f) == getfield(x.measurements[k], f)
            end
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    nothing
end

for file in readdir()
    if endswith(file, "jld") || endswith(file, "jld2") || endswith(file, ".confs")
        rm(file)
    end
end

@time @testset "DQMC" begin
    isfile("testfile.confs") && rm("testfile.confs")
    model = HubbardModel(4, 2, t = 1.7, U = 2.5)
    mc = DQMC(
        model, beta = 1.0, thermalization = 21, sweeps = 117, measure_rate = 1,
        recorder = BufferedConfigRecorder(DensityHirschField, "testfile.confs", rate = 1)
    )
    mc[:CDC] = charge_density_correlation(mc, model)
    run!(mc, verbose=false)

    save("testfile.jld", mc)
    x = load("testfile.jld")
    rm("testfile.jld")
    @test mc.field.conf == x.field.conf

    # Repeat these tests once with x being replayed rather than loaded
    test_dqmc(mc, x)

    # Check everything again with x being a replayed simulation
    x[:CDC] = charge_density_correlation(x, model)
    x.last_sweep = 0
    replay!(x, verbose=false)
    test_dqmc(mc, x)
    isfile("testfile.confs") && rm("testfile.confs")

    # Test resume

    # Run for 1s with known RNG
    Random.seed!(123)
    model = HubbardModel(2, 2, t = 1.7, U = 2.5)
    mc = DQMC(
        model, beta = 1.0, thermalization = 0, sweeps = 10_000_000, measure_rate = 100,
        recorder = BufferedConfigRecorder(DensityHirschField, "testfile.confs", rate = 100)
    )
    mc[:CDC] = charge_density_correlation(mc, model)

    state = run!(
        mc, verbose = false,
        safe_before = now() + Second(2),
        grace_period = Millisecond(0),
        resumable_filename = "resumable_testfile.jld2"
    )

    @test state == false
    cs = [mc.recorder[i] for i in 1:length(mc.recorder)]
    @assert length(cs) > 1 "No measurements have been taken. Test with more time!"
    L = length(cs)

    # Test whether safe file gets overwritten correctly
    mc, state = resume!(
        "resumable_testfile.jld2",
        verbose = false,
        safe_before = now() + Second(10),
        grace_period = Millisecond(0),
        overwrite = true,
        resumable_filename = "resumable_testfile.jld2"
    )

    @test state == false
    cs = [mc.recorder[i] for i in 1:length(mc.recorder)]
    @assert length(cs) - L > 1 "No new measurements have been taken. Test with more time!"
    @test isfile("resumable_testfile.jld2")
    isfile("testfile.confs") && rm("testfile.confs")

    # Test whether data from resumed simulation is correct
    Random.seed!(123)
    model = HubbardModel(2, 2, t = 1.7, U = 2.5)
    mc = DQMC(
        model, beta = 1.0, thermalization = 0, sweeps = 100length(cs), measure_rate = 100,
        recorder = BufferedConfigRecorder(DensityHirschField, "testfile.confs", rate = 100)
    )
    state = run!(mc, verbose = false)
    matches = true
    for i in 1:length(mc.recorder)
        matches = matches && (mc.recorder[i] == cs[i])
    end
    @test matches
    rm("resumable_testfile.jld2")
    isfile("testfile.confs") && rm("testfile.confs")
end


text:
using BenchmarkTools
using TextParse
VERSION >= v"0.7" && using Dates

const SUITE = BenchmarkGroup()

SUITE["util"] = BenchmarkGroup(["string", "unicode"])

our_lastindex(x) = VERSION >= v"0.7" ? lastindex(x) : endof(x)

float64str = "456.254"
float64strlen = our_lastindex(float64str)
negfloat64str = "-456.254"
negfloat64strlen = our_lastindex(negfloat64str)

intstr = "9823345"
intstrlen = our_lastindex(intstr)
negintstr = "-9823345"
negintstrlen = our_lastindex(negintstr)

SUITE["util"]["tryparsenext_base10_digit"] = BenchmarkGroup()
SUITE["util"]["tryparsenext_base10_digit"]["Float64"] = @benchmarkable TextParse.tryparsenext_base10_digit(Float64, $float64str, 3, $float64strlen)
SUITE["util"]["tryparsenext_base10_digit"]["Int64"] = @benchmarkable TextParse.tryparsenext_base10_digit(Int64, $intstr, 3, $intstrlen)

SUITE["util"]["tryparsenext_base10"] = BenchmarkGroup()
SUITE["util"]["tryparsenext_base10"]["Float64"] = @benchmarkable TextParse.tryparsenext_base10(Float64, $float64str, 3, $float64strlen)
SUITE["util"]["tryparsenext_base10"]["Int64"] = @benchmarkable TextParse.tryparsenext_base10(Int64, $intstr, 3, $intstrlen)

SUITE["util"]["tryparsenext_sign"] = BenchmarkGroup()
SUITE["util"]["tryparsenext_sign"]["nosign"] = BenchmarkGroup()
SUITE["util"]["tryparsenext_sign"]["nosign"]["Float64"] = @benchmarkable TextParse.tryparsenext_sign($float64str, 1, $float64strlen)
SUITE["util"]["tryparsenext_sign"]["nosign"]["Int64"] = @benchmarkable TextParse.tryparsenext_sign($intstr, 1, $intstrlen)
SUITE["util"]["tryparsenext_sign"]["neg"] = BenchmarkGroup()
SUITE["util"]["tryparsenext_sign"]["neg"]["Float64"] = @benchmarkable TextParse.tryparsenext_sign($negfloat64str, 1, $negfloat64strlen)
SUITE["util"]["tryparsenext_sign"]["neg"]["Int64"] = @benchmarkable TextParse.tryparsenext_sign($negintstr, 1, $negintstrlen)

whitespacestring = "abc  de"
nowhitespacestring = "abcde"

SUITE["util"]["eatwhitespaces"] = BenchmarkGroup()
SUITE["util"]["eatwhitespaces"]["withwhitespace"] = @benchmarkable TextParse.eatwhitespaces($whitespacestring, 4, $(lastindex(whitespacestring)))
SUITE["util"]["eatwhitespaces"]["nowhitespacestring"] = @benchmarkable TextParse.eatwhitespaces($whitespacestring, 4, $(lastindex(whitespacestring)))

newlinestring = "ab\r\n\r\r"

SUITE["util"]["eatnewlines"] = BenchmarkGroup()
SUITE["util"]["eatnewlines"]["default"] = @benchmarkable TextParse.eatnewlines($newlinestring, 3, $(lastindex(newlinestring)))

SUITE["util"]["getlineend"] = BenchmarkGroup()
SUITE["util"]["getlineend"]["default"] = @benchmarkable TextParse.getlineend($newlinestring)

percentagestring = "35.35%"
percentagestringlen = our_lastindex(percentagestring)
somestring = "foo something,"
somestringlen = our_lastindex(somestring)
somequotedstring =  "\"Owner 2 ”Vicepresident\"\"\""
somequotedstringlen = our_lastindex(somequotedstring)

longfloat64str = "2344345.1232353459389238738435"
longfloat64strlen = our_lastindex(longfloat64str)

tok = TextParse.DateTimeToken(DateTime, dateformat"yyyy-mm-dd HH:MM:SS")
opts = TextParse.LocalOpts('y', false, '"', '\\', false, false)
datetimestr = "1970-02-02 02:20:20"
datetimestrlen = our_lastindex(datetimestr)

SUITE["util"]["tryparsenext"] = BenchmarkGroup()
SUITE["util"]["tryparsenext"]["NumericFloat64"] = @benchmarkable TextParse.tryparsenext($(TextParse.Numeric(Float64)), $float64str,1,$float64strlen)
SUITE["util"]["tryparsenext"]["LongNumericFloat64"] = @benchmarkable TextParse.tryparsenext($(TextParse.Numeric(Float64)), $longfloat64str,1,$longfloat64strlen)
SUITE["util"]["tryparsenext"]["UInt64"] = @benchmarkable TextParse.tryparsenext($(TextParse.Numeric(UInt64)), $intstr,1,$intstrlen)
SUITE["util"]["tryparsenext"]["NegInt64"] = @benchmarkable TextParse.tryparsenext($(TextParse.Numeric(Int64)), $negintstr,1,$negintstrlen)
SUITE["util"]["tryparsenext"]["Percentage"] = @benchmarkable TextParse.tryparsenext($(TextParse.Percentage()), $percentagestring,1,$percentagestringlen, TextParse.default_opts)
SUITE["util"]["tryparsenext"]["StringToken"] = @benchmarkable TextParse.tryparsenext($(TextParse.StringToken(String)), $somestring,1,$somestringlen, TextParse.default_opts)
SUITE["util"]["tryparsenext"]["DateTimeToken"] = @benchmarkable TextParse.tryparsenext($tok, $datetimestr,1,$datetimestrlen, $opts)
SUITE["util"]["tryparsenext"]["QuotedStringToken"] = @benchmarkable TextParse.tryparsenext($(Quoted(String,quotechar='"', escapechar='"')), $somequotedstring)

somefieldstring = " 12,3"
f = TextParse.fromtype(Int)
SUITE["util"]["tryparsenext"]["Field"] = @benchmarkable TextParse.tryparsenext($(TextParse.Field(f)), $somefieldstring)


text:
using Prototypes
using Test, BenchmarkTools

struct A <: Prototyped
    values::Dict{Symbol, Any}
    __proto
    A(prototype=Object()) = new(Dict(), prototype)
end

struct B{TProto} <: Prototyped
    values::Dict{Symbol, Any}
    __proto::TProto
    B(prototype=Object()) = new{typeof(prototype)}(Dict(), prototype)
end

@testset "basics" begin
    a1 = A()
    a1.vala = 42
    @test_throws ErrorException a1.valx
    @test a1.vala === 42
    a2 = A(a1)
    a2.vala = 1
    @test a2.vala === 1
    b1 = B(a1)
    b1.valb = "b1"
    b1.valb1 = "b11"
    b2 = B(b1)
    b2.valb = "b2"
    @test b1.valb === "b1"
    @test b1.valb1 === "b11"
    @test b1.vala === 42
    @test_throws ErrorException b1.valx
    @test b2.valb === "b2"
    @test b2.valb1 === "b11"
    @test b2.vala === 42
    @test_throws ErrorException b2.valx

    a1.valx = 20
    @test b2.valx === 20
    @test b1.valx === 20
    @test a2.valx === 20
    @test a1.valx === 20
    b1.valx = 30
    @test b2.valx === 30
    @test b1.valx === 30
    @test a2.valx === 20
    @test a1.valx === 20
end

@testset "benchmarks" begin
    a1 = A()
    a1.vala = 42
    b1 = B(a1)
    b2 = B(b1)
    b3 = B(b2)
    b4 = B(b3)

    @btime $a1.vala
    @btime $b1.vala
    @btime $b2.vala
    @btime $b3.vala
    @btime $b4.vala

    a1.a1 = 1
    b1.b1 = 1
    b2.b2 = 1
    b3.b3 = 1
    b4.b4 = 1
    @test (@btime $a1.vala) === 42
    @btime $b1.vala
    @btime $b2.vala
    @btime $b3.vala
    @btime $b4.vala
end


text:
# Bayesian Logistic Regression using the Polya-Gamma augmentation.
#

using LinearAlgebra

export BinaryLogisticRegression
export LogisticRegression
export predict

#######################################################################
# Binary logistic regression
#
# NOTE:
# We discourage the use of BinaryLogisticRegression model directly.
# Instead, you should use the general LogisticRegression model (which
# uses the BinaryLogisticRegression itself).

struct BinaryLogisticRegression <: Model
    β::ConjugateParameter{<:Normal}
    hasbias::Bool

    function BinaryLogisticRegression(β::ConjugateParameter{<:Normal}, hasbias::Bool)
        model = new(β, hasbias)
        β.stats = data -> stats_β(model, data...)
        return model
    end
end

function BinaryLogisticRegression(μ₀::Vector{T}, Σ₀::Matrix{T},
                                  μ::Vector{T}, Σ::Matrix{T},
                                  hasbias::Bool) where T <: AbstractFloat
    BinaryLogisticRegression(ConjugateParameter(Normal(μ₀, Σ₀), Normal(μ, Σ)),
                             hasbias)
end

function BinaryLogisticRegression(T; inputdim::Integer, initσ::Real = 0.0,
                                  pseudocounts::Real = 1.0, hasbias::Bool = true)
    dim = hasbias ? inputdim + 1 : inputdim

    # Prior parameters
    μ₀ = zeros(T, dim)
    Σ₀ = T(1/pseudocounts) * Matrix{T}(I, dim, dim)

    # Initial paremeterization of the posterior
    μ = μ₀ .+ T(initσ) .* randn(T)
    Σ = deepcopy(Σ₀)

    BinaryLogisticRegression(μ₀, Σ₀, μ, Σ, hasbias)
end

function BinaryLogisticRegression(;inputdim::Integer, initσ::Real = 0.0,
                                  pseudocounts::Real = 1.0, hasbias::Bool = true)
    BinaryLogisticRegression(Float64, inputdim = inputdim, initσ = initσ,
                             pseudocounts = pseudocounts, hasbias = hasbias)
end

function regressors(model::BinaryLogisticRegression, X::Matrix{T}) where T <: AbstractFloat
    if model.hasbias
        # We add an extra dimension always set to "1" to incorporate a bias
        # parameter to the model.
        return vcat(X, ones(T, 1, size(X, 2)))
    end
    return X
end

# Return the Polya-Gamma posterior of the augmented
# logistic function
function _augment_posterior(ψ²::Vector{T}) where T <: AbstractFloat
    b = ones(T, length(ψ²))
    qω = PolyaGamma(b)
    update!(qω, -.5 * ψ²)
    return qω
end

# Return the sufficient statistics of the log likelihood and the
# the expected statistics of β.
function _precompute(model::BinaryLogisticRegression, X::Matrix{T1},
                     z::Vector{T2},
                     Nₖ::Union{Vector{T2}, Nothing} = nothing) where {T1 <: AbstractFloat,
                                                                      T2 <: Real}
    X̂ = regressors(model, X)

    # Quadratic expansion of the regressors
    D, N = size(X̂)
    X̂X̂ᵀ = reshape(X̂, D, 1, N) .* reshape(X̂, 1, D, N)
    vec_X̂X̂ᵀ = reshape(X̂X̂ᵀ, :, N)

    # Expectation of the regressors: E[ T(β) ] = E[ ( β, vec(ββᵀ) )]
    E_Tβ = gradlognorm(model.β.posterior)
    vec_E_ββᵀ = E_Tβ[D+1:end]

    # Augmentation of the logistic function
    pω = PolyaGamma{T1, size(X̂, 2)}()
    ψ² = vec_X̂X̂ᵀ' * vec_E_ββᵀ
    qω = _augment_posterior(ψ²)
    E_ω = mean(qω)

    zstats = isnothing(Nₖ) ? (z .- 0.5) : z .- Nₖ * 0.5

    stats = vcat(reshape(zstats, 1, N) .* X̂,
                 -.5 .* reshape(E_ω, 1, :) .* vec_X̂X̂ᵀ)

    stats, E_Tβ, pω, qω
end

function (model::BinaryLogisticRegression)(X::Matrix{T1},
                                           z::Vector{T2},
                                           Nₖ::Union{Vector{T2}, Nothing} = nothing) where {T1 <: AbstractFloat,
                                                                                T2 <: Real}
    stats, E_Tβ, pω, qω = _precompute(model, X, z, Nₖ)

    # KL divegerence betwen qω and pω for each sample.
    q_η, p_η = naturalparam(qω), naturalparam(pω)
    KL = lognorm(pω, perdim = true) - lognorm(qω, perdim = true) - (p_η .- q_η) .* gradlognorm(qω)

    stats' * E_Tβ .- KL .- log(2)
end

function stats_β(model::BinaryLogisticRegression,
                       X::Matrix{T1},
                       z::Vector{T2},
                       Nₖ::Union{Vector{T2}, Nothing} = nothing) where {T1 <: AbstractFloat,
                                                                        T2 <: Real}
    stats, E_Tβ, _, _ = _precompute(model, X, z, Nₖ)
    stats
end

function predict(model::BinaryLogisticRegression,
                 X::Matrix{T}; a::T = 0.368) where T <: AbstractFloat
    # Following https://arxiv.org/pdf/1703.00091.pdf, we use the
    # approximation:
    # ⟨σ(ψ)⟩ ≈ σ(μ / √(1 + a * σ²))
    # where a = 0.368
    # In Bishop 4.5.2 a very similar approximation uses a = π / 8 ≈ 0.392...

    μ, Σ = stdparam(model.β.posterior)
    X̂ = regressors(model, X)

    σ²s = [sum((Σ .* x̂') .* x̂) for x̂ in eachcol(X̂)]
    ψs = (X̂' * μ) ./ sqrt.(1 .+ a * σ²s)
    1 ./ (1 .+ exp.(-ψs))
end

#######################################################################
# K-class Logistic Regression (LR) based on the stick-breaking process
#
# References
# ----------
# [1] Dependent Multinomial Modes Made Easy: Stick-Breaing with the
#     Polya-Gamma Augmentation (https://arxiv.org/pdf/1506.05843.pdf)

struct LogisticRegression{K} <: Model
    # The LR is defined by K - 1 binary LR models (see eq. (7) in [1])
    stickbreaking::Array{BinaryLogisticRegression}

    function LogisticRegression(sb::Array{BinaryLogisticRegression})
        model = new{length(sb) + 1}(sb)

        # We override the stats
        for (k, binarylr) in enumerate(sb)
            binarylr.β.stats = data -> stats_β(model, k, data...)
        end

        return model
    end
end

function LogisticRegression(μ₀::Vector{T}, Σ₀::Matrix{T}, μ::Vector{T},
                            Σ::Matrix{T}, hasbias::Bool; nclasses::Integer) where T <: AbstractFloat
    LogisticRegression([BinaryLogisticRegression(μ₀, Σ₀, deepcopy(μ), deepcopy(Σ), hasbias)
                        for i in 1:nclasses-1])
end

function LogisticRegression(T; inputdim::Integer, nclasses::Integer,
                            initσ::Real = 0.0, pseudocounts::Real = 1.0,
                            hasbias::Bool = true)
    dim = hasbias ? inputdim + 1 : inputdim

    # Prior parameters
    μ₀ = zeros(T, dim)
    Σ₀ = T(1/pseudocounts) * Matrix{T}(I, dim, dim)

    # Initial paremeterization of the posterior
    μ = μ₀ .+ T(initσ) .* randn(T)
    Σ = deepcopy(Σ₀)

    LogisticRegression(μ₀, Σ₀, μ, Σ, hasbias, nclasses = nclasses)
end

function LogisticRegression(;inputdim::Integer, nclasses::Integer,
                            initσ::Real = 0.0, pseudocounts::Real = 1.0,
                            hasbias::Bool = true)
    LogisticRegression(Float64, inputdim = inputdim, nclasses = nclasses,
                       initσ = initσ, pseudocounts = pseudocounts,
                       hasbias = hasbias)
end

getconjugateparams(model::LogisticRegression) = [m.β for m in model.stickbreaking]

function _encode(::LogisticRegression{K}, z::Vector{T}) where {K, T <: Integer}
    onehot = zeros(T, K - 1, length(z))
    for i = 1:K-1
        onehot[i, z .== i] .= 1
    end

    # See eq. (6) in [1]. In our case, N = 1
    Nₖ = zeros(T, K - 1, length(z))
    for i = 1:K-1
        Nₖ[i, z .>= i] .= 1
    end

   onehot, Nₖ
end

function (model::LogisticRegression)(X::Matrix{T1}, z::Vector{T2},) where {T1 <: AbstractFloat,
                                                                           T2 <: Integer}
    onehot, Nₖ = _encode(model, z)
    retval = zeros(T1, size(X, 2))
    for (k, m) in enumerate(model.stickbreaking)
        retval[z .>= k] .+= m(X[:, z .>= k], onehot[k, z .>= k], Nₖ[k, z .>= k])
    end
    retval
end

function stats_β(model::LogisticRegression, k::Integer,
                         X::Matrix{T1},
                         z::Vector{T2}) where {T1 <: AbstractFloat, T2 <: Real}
    onehot, Nₖ = _encode(model, z)
    stats_β(model.stickbreaking[k], X[:, z .>= k], onehot[k, z .>= k],
                  Nₖ[k, z .>= k])
end

function predict(model::LogisticRegression{K}, X::Matrix{T}) where {K, T <: AbstractFloat}
    retval = zeros(T, K, size(X, 2))
    for k in 1:K
        residual = dropdims(sum(retval[1:k, :], dims = 1), dims = 1)
        if k < K
            retval[k, :] = predict(model.stickbreaking[k], X) .* (1 .- residual)
        else
            retval[k, :] = (1 .- residual)
        end
    end
    retval
end

max_stars_repo_name: fnlandini/BayesianModels

max_issues_repo_name: fnlandini/BayesianModels

max_forks_repo_name: fnlandini/BayesianModels


text:
save(list = ls(all.names = TRUE), file = "/media/u01/analytics/scoring/CDW5_792/RConvert/renv.RData", envir = .GlobalEnv)

# =========== JULIA ===========

    R""" load("/media/u01/analytics/scoring/CDW5_792/RConvert/renv.RData")
         ls()
     """

# ------- CFG -------
#function fcmp(l1:)
#end
chk=OrderedDict()
        chk[:P2_Competitor] = cfg[:P2_Competitor] == rcopy("P2_Competitor")
        chk[:pvalue_lvl] = cfg[:pvalue_lvl] == rcopy("pvalue_lvl")
        chk[:exposed_flag_var] = string(cfg[:exposed_flag_var]) == rcopy("exposed_flag_var")
        chk[:random_demos] = sort(ASCIIString[string(x) for x in cfg[:random_demos]]) == sort(rcopy("random_demos"))
        chk[:random_campaigns] = sort(ASCIIString[string(x) for x in cfg[:random_campaigns]]) == sort(rcopy("random_campaigns"))
        #sort(ASCIIString[string(x) for x in cfg[:random_campaigns]]) == sort(rcopy("dropvars"))
        if (typeof(rcopy("dropvars")) == ASCIIString) & (length(cfg[:dropvars]) == 1)
            chk[:dropvars] = string(cfg[:dropvars][1]) == rcopy("dropvars")
        else
            chk[:dropvars] = sort(ASCIIString[string(x) for x in cfg[:dropvars]]) == sort(rcopy("dropvars"))
        end
        chk[:scoring_vars] = sort(ASCIIString[string(x) for x in cfg[:scoring_vars]]) == sort(rcopy("scoring_vars"))
        chk[:all_mandatory_vars] = sort(ASCIIString[string(x) for x in cfg[:all_mandatory_vars]]) == sort(rcopy("all_mandatory_vars"))
        chk[:allrandoms] = sort(ASCIIString[string(x) for x in cfg[:allrandoms]]) == sort(rcopy("allrandoms"))
        chk[:ProScore] = string(cfg[:ProScore]) == rcopy("ProScore")
        if rcopy("exists('random_demos_factor')")
            chk[:random_demos_factor] = sort(ASCIIString[string(x) for x in cfg[:random_demos_factor]]) == sort(rcopy("random_demos_factor"))
        else
            chk[:random_demos_factor] =  "random_demos_factor doesn't exist"
        end
        chk[:num_products] = cfg[:num_products] == rcopy("num_products")
        #addedvars=[:isO,:whyO,:data_NB_NE_B,:data_B_E_NB,:pen_reduction,:occ_reduction,:dolocc_reduction]
        chk[:xVarsDemos] =  sort(ASCIIString[string(x) for x in cfg[:xVarsDemos]]) == sort(rcopy("xVarsDemos"))
        #setdiff(sort(rcopy("xVarsDemos")),sort(ASCIIString[string(x) for x in cfg[:xVarsDemos]]))
        chk[:xVarsPost] =  sort(ASCIIString[string(x) for x in cfg[:xVarsPost]]) == sort(rcopy("initial_vars_to_exclude"))
        #setdiff(sort(rcopy("initial_vars_to_exclude")),sort(ASCIIString[string(x) for x in cfg[:xVarsPost]]))
        chk[:iVarsPREPOS] = sort(ASCIIString[string(x) for x in cfg[:iVarsPREPOS]]) == sort(rcopy("PREPOS_vars_to_include"))

        chk[:xVarsP0] =  sort(ASCIIString[string(x) for x in cfg[:xVarsP0]]) == sort(rcopy("P0_vars_to_exclude"))
        #setdiff(sort(rcopy("P0_vars_to_exclude")),sort(ASCIIString[string(x) for x in cfg[:xVarsP0]]))
        #setdiff(sort(ASCIIString[string(x) for x in cfg[:xVarsP0]]), sort(rcopy("P0_vars_to_exclude")))

        chk[:xVarsReports] = sort(ASCIIString[string(x) for x in cfg[:xVarsReports]]) == sort(rcopy("vars_to_exclude"))
        #setdiff(sort(rcopy("vars_to_exclude")),sort(ASCIIString[string(x) for x in cfg[:xVarsReports]]))
        #setdiff(sort(ASCIIString[string(x) for x in cfg[:xVarsReports]]),sort(rcopy("vars_to_exclude")))

        #chk[:xvars] = sort(ASCIIString[string(x) for x in cfg[:xvars]]) == sort(rcopy("ALL_vars_to_exclude"))
        ##setdiff(sort(rcopy("ALL_vars_to_exclude")),sort(ASCIIString[string(x) for x in cfg[:xvars]]))
        ##setdiff(sort(ASCIIString[string(x) for x in cfg[:xvars]]),sort(rcopy("ALL_vars_to_exclude")))
        #:xvars
        #:ivars
        #chk[:ALL_vars_to_exclude] = sort(ASCIIString[string(x) for x in cfg[:ALL_vars_to_exclude]]) == sort(rcopy("ALL_vars_to_exclude"))
        ##setdiff(sort(rcopy("ALL_vars_to_exclude")),sort(ASCIIString[string(x) for x in cfg[:ALL_vars_to_exclude]]))
        ##setdiff(sort(ASCIIString[string(x) for x in cfg[:ALL_vars_to_exclude]]),sort(rcopy("ALL_vars_to_exclude")))

        chk[:negativevars] =  sort(ASCIIString[string(x) for x in cfg[:negativevars]]) == sort(rcopy("negativevars"))
        #setdiff(sort(rcopy("negativevars")),sort(ASCIIString[string(x) for x in cfg[:negativevars]]))
        #setdiff(sort(ASCIIString[string(x) for x in cfg[:negativevars]]), sort(rcopy("negativevars")))

        chk[:positivevars] =  sort(ASCIIString[string(x) for x in cfg[:positivevars]]) == sort(rcopy("positivevars"))
        #setdiff(sort(rcopy("positivevars")),sort(ASCIIString[string(x) for x in cfg[:positivevars]]))
        #setdiff(sort(ASCIIString[string(x) for x in cfg[:positivevars]]), sort(rcopy("positivevars")))
        #:ALL_vars_to_exclude
        #:negativevars
        #:positivevars
        initial_data=rcopy("initial_data")
        chk[:features] = sort(names(df_data)) == sort(names(initial_data))
        #setdiff(sort(names(initial_data)),sort(names(df_data)))
        #setdiff(sort(names(df_data)) , sort(names(initial_data))  )

# ---- DF_DATA Check -----
save(list = ls(all.names = TRUE), file = "/media/u01/analytics/scoring/CDW5_792/RConvert/renv.RData", envir = .GlobalEnv)

R""" load("/media/u01/analytics/scoring/CDW5_792/RConvert/renv.RData")
    ls()
 """
greparr=names(df_in)
function grep(v::AbstractString,xvars::Array{Symbol}=greparr)   #function grep(v::AbstractString) filter(x->ismatch(Regex("$(v)"), string(x)), vars) end
     filter(x->ismatch(Regex("$(v)"), string(x)), xvars)
end
function grep(vi::Symbol,xvars::Array{Symbol}=greparr)   #function grep(v::AbstractString) filter(x->ismatch(Regex("$(v)"), string(x)), vars) end
    v = string(vi)
    filter(x->ismatch(Regex("$(v)"), string(x)), xvars)
end
finaldata = rcopy("finaldata")
initial_data = finaldata
initial_data = rcopy("initial_data")

# ------------------------------------------------------------------------------------------------------
#@rget finaldata;
jlst = sort(setdiff(names(df_data),[:data_B_E_NB,:data_NB_NE_B,:occ_reduction,:dolocc_reduction ,:pen_reduction,]))

setdiff(names(df_data),names(finaldata))
setdiff(names(df_data),names(initial_data))
#setdiff(names(df_in),names(rcopy("finaldata")))
setdiff(names(df_in),names(finaldata))
grep(:exposed_flag,names(df_in))
grep(:exposed_flag,names(finaldata))

println("Only In df_in:::\n",setdiff(names(df_in),names(finaldata))  )
println("Only In finaldata:::\n",setdiff(names(finaldata), names(df_in)))

# --- Data ----
println("df_data : ",nrow(df_data),"  ---->  initial_data : ",nrow(initial_data))
#df_data[df_data[:number_of_children_in_living_Un].>=4,:number_of_children_in_living_Un] = 4  # aggregate #of children for 4+

#NOTTE: remove observations;
#            291 - remove person_1_gender
#            298 - remove HH exposed/non-exposed
#            look for dropleveles - not sureif has impact
#df_data=df_in
# :number_of_children_in_living_Un  # 4+ values
#  df_data[df_data[:number_of_children_in_living_Un].>=4,:number_of_children_in_living_Un] = 4
function dict2str(di::Dict)
    od=Dict()
    for (k, v) in di
        od[string(k)] = string(v)
    end
    return od
end
function cm()
    di=dict2str(countmap(df_in[:number_of_children_in_living_Un]))
    dd=dict2str(countmap(df_data[:number_of_children_in_living_Un]))
    fd=dict2str(countmap(finaldata[:number_of_children_in_living_Un]))
    id=dict2str(countmap(initial_data[:number_of_children_in_living_Un]))
    #for k in sort(unique([replace(string(x),"+","") for x in union(keys(dd),keys(fd),keys(id))]))

    dfo=DataFrame(key=AbstractString[],df_in=AbstractString[],df_data=AbstractString[],finaldata=AbstractString[],initial_data=AbstractString[])
    for k in sort(unique([string(x) for x in union(keys(di),keys(dd),keys(fd),keys(id))]))
        iv=AbstractString[k,get(di,k,""),get(dd,k,""),get(fd,k,""),get(id,k,"")]
        push!(dfo,iv)
        #push!(df, [3  6])
        #push!(dfo,[get(dd,k,""),get(fd,k,""),get(id,k,"")])
        #println(":::",iv)
    end
    return dfo
end
cm()

#   :group   # values and counts/proportions
function cm()
    di=dict2str(countmap(df_in[ df_data[:isO].==false, :group]))
    fd=dict2str( countmap(finaldata[:group]) )
    dfo=DataFrame(key=AbstractString[],df_in=AbstractString[],finaldata=AbstractString[])
    for k in sort(unique([string(x) for x in union(keys(di),keys(fd),)]))
        iv=AbstractString[k,get(di,k,""),get(fd,k,"")]
        push!(dfo,iv)
        #push!(df, [3  6])
        #push!(dfo,[get(dd,k,""),get(fd,k,""),get(id,k,"")])
        #println(":::",iv)
    end
    return dfo
end
cm()

# L_290 remove HHs with no gender info  #finaldata <- finaldata[!person_1_gender=="U",]
function cm()
    di=dict2str(countmap(df_in[  df_data[:isO].==false, :person_1_gender]))
    fd=dict2str( countmap(finaldata[:person_1_gender]) )
    dfo=DataFrame(key=AbstractString[],df_in=AbstractString[],finaldata=AbstractString[])
    for k in sort(unique([string(x) for x in union(keys(di),keys(fd),)]))
        iv=AbstractString[k,get(di,k,""),get(fd,k,"")]
        push!(dfo,iv)
        #push!(df, [3  6])
        #push!(dfo,[get(dd,k,""),get(fd,k,""),get(id,k,"")])
        #println(":::",iv)
    end
    return dfo
end
cm()

# # aggregate U and L levels of hh income  #finaldata <- finaldata[estimated_hh_income=="U",estimated_hh_income:="L"]
function cm()
    di=dict2str(countmap(df_in[  df_data[:isO].==false, :estimated_hh_income]))
    fd=dict2str( countmap(finaldata[:estimated_hh_income]) )
    dfo=DataFrame(key=AbstractString[],df_in=AbstractString[],finaldata=AbstractString[])
    for k in sort(unique([string(x) for x in union(keys(di),keys(fd),)]))
        iv=AbstractString[k,get(di,k,""),get(fd,k,"")]
        push!(dfo,iv)
        #push!(df, [3  6])
        #push!(dfo,[get(dd,k,""),get(fd,k,""),get(id,k,"")])
        #println(":::",iv)
    end
    return dfo
end
cm()

function cm()
    o=Dict()
    for r in cfg[:random_campaigns]
        #di = countmap(df_data[df_data[:isO].==false, r])
        di = dict2str(countmap(df_data[df_data[:isO].==false, r]) )
        fd = dict2str( countmap(finaldata[r]) )
        dfo=DataFrame(key=AbstractString[],df_data=AbstractString[],finaldata=AbstractString[])
        for k in sort(unique([string(x) for x in union(keys(di),keys(fd),)]))
            iv=AbstractString[k,get(di,k,""),get(fd,k,"")]
            push!(dfo,iv)
        end
        o[r] = dfo
    end
    return o
end
o=cm()

nrow(df_data[df_data[:isO].==false, : ])
nrow(initial_data)

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
function create_ols_approximation(y::Array{<:Real,1}, x::Array{<:Real,1}, base_x::Real = 0.0, degree::Integer = 1, intercept::Bool = true)
    obs = length(y)
    if degree < 0
        error("Cannot approximate with OLS with a degree that is negative")
    end
    x = x .- base_x
    if intercept
        X = ones(obs)
        for i in 1:degree
            X = hcat(X, (x .^ i))
        end
    else
        X = x
        for i in 2:degree
            X = hcat(X, (x .^ i))
        end
    end

    lm1 = fit(LinearModel,  hcat(X), y)
    beta = lm1.pp.beta0
    func_array = Array{PE_Function,1}(undef,convert(Int, intercept) + degree)
    if intercept
        func_array[1] = PE_Function(beta[1], 0.0, base_x, 0)
    end
    for d in 1:degree
        func_array[d+Integer(intercept)] = PE_Function(beta[d+Integer(intercept)], 0.0, base_x, d)
    end
    return Sum_Of_Functions(func_array)
end

function create_ols_approximation(y::Array{<:Real,1}, x::Union{Array{DateTime,1},Array{Date,1},Array{Union{Date,DateTime},1}}, base_x::Union{Date,DateTime} = global_base_date, degree::Integer = 1, intercept::Bool = true)
    base = years_from_global_base.(base_x)
    xx   = years_from_global_base.(x)
    return create_ols_approximation(y, xx, base, degree, intercept)
end

function get_cholesky_coefficients(chebyshev::Sum_Of_Functions, y::Array{<:Real,1}, unnormalised_nodes::Array{<:Real,1})
    chebyshev_on_nodes = evaluate.(Ref(chebyshev), unnormalised_nodes)
    a = sum(y .* chebyshev_on_nodes) / sum(chebyshev_on_nodes .^ 2)
    return a
end

function create_chebyshev_approximation(func::Function, nodes::Integer, degree::Integer, left::Real, right::Real)
    # This is all after Algorithm 6.2 from Judd (1998) Numerical Methods in Economics.
    if nodes <= degree
        error("Need to have more nodes than degree to use a chebyshev approximation")
    end
    k = 1:nodes
    unnormalised_nodes = -cos.( (((2 .* k) .- 1) ./ (2 * nodes)) .* pi    )
    normalised_nodes = ((unnormalised_nodes .+ 1) .* ((right-left)/2)) .+ left
    y = func.(normalised_nodes)
    chebyshevs = get_chevyshevs_up_to(degree, true)
    a = get_cholesky_coefficients.(chebyshevs, Ref(y), Ref(unnormalised_nodes))
    transformed_chebyshevs = convert_to_linearly_rescale_inputs.(chebyshevs, (right-left)/2, +(right+left)/2)
    # Note that these alpha and beta parameters in the convert_to_linearly_rescale_inputs function differ from those in Judd because those did not work in this context.
    all_terms = a .* transformed_chebyshevs
    final_func = Sum_Of_Functions(all_terms)
    return final_func
end