text:
module ScientificTypes

# Dependencies
using Reexport
@reexport using ScientificTypesBase
using Tables
using CategoricalArrays
using ColorTypes
using PrettyTables
using Dates
import Distributions

import StatisticalTraits: info

# re-export from StatisticalTraits
export info

# exports
export coerce, coerce!, autotype, schema

# -------------------------------------------------------------
# Abbreviations

const ST   = ScientificTypesBase
const Arr  = AbstractArray
const CArr = CategoricalArray
const Cat  = CategoricalValue

# Indicate the convention, see init.jl where it is set.
struct DefaultConvention <: Convention end

include("init.jl")

# -------------------------------------------------------------
# Includes

include("coerce.jl")
include("schema.jl")
include("autotype.jl")

include("convention/scitype.jl")
include("convention/coerce.jl")

end # module


text:
function read_request(channel::gRPCChannel, controller::gRPCController, services)
    connection = channel.session
    evt = Session.take_evt!(connection)

    if isa(evt, EvtGoaway)
        debug_log(controller, "received EvtGoaway")
        close(channel)
        return (nothing, nothing, nothing)
    end

    if !isa(evt, EvtRecvHeaders)
        debug_log(controller, "unexpected event $evt")
        close(channel)
        return (nothing, nothing, nothing)
    end

    channel.stream_id = evt.stream_identifier

    headers = evt.headers
    method = headers[":method"]
    path = headers[":path"]
    pathcomps = split(path, "/"; keep=false)
    if length(pathcomps) != 2
        debug_log(controller, "unexpected path $path")
        close(channel)
        return (nothing, nothing, nothing)
    end

    servicename, methodname = split(path, "/"; keep=false)

    if evt.is_end_stream
        data = UInt8[]
    else
        data_evt = Session.take_evt!(connection)
        data = data_evt.data
    end

    debug_log(controller, "request: $method $path")
    debug_log(controller, "stream: $(channel.stream_id)")
    debug_log(controller, "service: $servicename")
    debug_log(controller, "method: $methodname")
    debug_log(controller, "data: $(length(data)) bytes")

    service = services[servicename]
    method = find_method(service, methodname)
    request_type = get_request_type(service, method)
    request = request_type()
    from_delimited_message_bytes(data, request)

    service, method, request
end

function write_response(channel::gRPCChannel, controller::gRPCController, response)
    sending_headers = Headers(":status" => "200")
    data_buff = to_delimited_message_bytes(response)

    Session.put_act!(channel.session, Session.ActSendHeaders(channel.stream_id, sending_headers, false))
    Session.put_act!(channel.session, Session.ActSendData(channel.stream_id, data_buff, true))
    nothing
end

# gRPC server implementation
#
type gRPCServer
    sock::TCPServer
    services::Dict{String, ProtoService}
    run::Bool
    debug::Bool

    gRPCServer(services::Tuple{ProtoService}, ip::IPv4, port::Integer) = gRPCServer(services, listen(ip, port))
    gRPCServer(services::Tuple{ProtoService}, port::Integer) = gRPCServer(services, listen(port))
    function gRPCServer(services::Tuple{ProtoService}, sock::TCPServer)
        svcdict = Dict{String,ProtoService}()
        for svc in services
            svcdict[svc.desc.name] = svc
        end
        new(sock, svcdict, true, false)
    end
end

# TODO: close all channels, wait for/interrupt processors
close(srvr::gRPCServer) = close(srvr.sock)

function process(controller::gRPCController, srvr::gRPCServer, channel::gRPCChannel)
    debug_log(controller, "start processing channel")
    try
        while(!channel.session.closed)
            service, method, request = read_request(channel, controller, srvr.services)
            (service === nothing) && continue

            response = call_method(service, method, controller, request)
            #debug_log(controller, "response: $response")
            write_response(channel, controller, response)
        end
    catch ex
        debug_log(controller, "channel stopped with exception $ex")
    end
    # TODO: close channel if not closed, remove reference from srvr
    debug_log(controller, "stopped processing channel")
end

function run(srvr::gRPCServer)
    controller = gRPCController(srvr.debug)
    try
        while(srvr.run)
            buffer = accept(srvr.sock)
            connection = Session.new_connection(buffer; isclient=false)

            # TODO: keep channel reference in server
            channel = gRPCChannel(connection)
            @async process(controller, srvr, channel)
        end
    catch ex
        debug_log(controller, "server stopped with exception $ex")
    end
    debug_log(controller, "stopped server")
end

ext: jl


text:
@testset "Principal Value (MRP)" begin
    for i in 1:1000
        p = MRP(5.0 * randn(), 5.0 * randn(), 5.0 * randn())
        p_prin = principal_value(p)
        @test p_prin ≈ p
        @test (p_prin.x^2 + p_prin.y^2 + p_prin.z^2) ≤ 1
    end
end

@testset "Principal Value (QuatRotation)" begin
    for i in 1:1000
        qq = rand(QuatRotation)
        qq_prin = principal_value(qq)
        @test 0.0 ≤ real(qq_prin.q)
        @test qq_prin ≈ qq
    end
end

@testset "Principal Value (Angle Axis)" begin
    for i in 1:1000
        aa = AngleAxis(100.0 * randn(), randn(), randn(), randn())
        aa_prin = principal_value(aa)
        @test 0.0 ≤ aa_prin.theta
        @test aa_prin ≈ aa
    end
end

@testset "Principal Value (Rotation Vector)" begin
    for i in 1:1000
        rv = RotationVec(100.0 * randn(), 100.0 * randn(), 100.0 * randn())
        rv_prin = principal_value(rv)
        @test rotation_angle(rv_prin) ≤ π
        @test rv_prin ≈ rv
    end
    rv = RotationVec(0.0, 0.0, 0.0)
    rv_prin = principal_value(rv)
    @test rotation_angle(rv_prin) ≤ π
    @test rv_prin ≈ rv
end

@testset "Principal Value (Rodrigues Parameters)" begin
    for i in 1:1000
        rv = RodriguesParam(100.0 * randn(), 100.0 * randn(), 100.0 * randn())
        rv_prin = principal_value(rv)
        @test rv_prin ≈ rv
        @test Rotations.params(rv_prin) ≈ Rotations.params(rv)
    end
    rv = RodriguesParam(0.0, 0.0, 0.0)
    rv_prin = principal_value(rv)
    @test rv_prin ≈ rv
    @test Rotations.params(rv_prin) ≈ Rotations.params(rv)
end

@testset "Principal Value ($(rot_type))" for rot_type in (RotX, RotY, RotZ, Angle2d)
    for i in 1:1000
        r = rot_type(100.0 * randn())
        r_prin = principal_value(r)
        @test -π ≤ r_prin.theta ≤ π
        @test r_prin ≈ r
    end
end

@testset "Principal Value ($(rot_type))" for rot_type in (RotXY, RotYX, RotZX, RotXZ, RotYZ, RotZY)
    for i in 1:1000
        r = rot_type(100.0 * randn(), 100.0 * randn())
        r_prin = principal_value(r)
        @test -π ≤ r_prin.theta1 ≤ π
        @test -π ≤ r_prin.theta2 ≤ π
        @test r_prin ≈ r
    end
end

@testset "Principal Value ($(rot_type))" for rot_type in (RotXYX, RotYXY, RotZXZ, RotXZX, RotYZY, RotZYZ, RotXYZ, RotYXZ, RotZXY, RotXZY, RotYZX, RotZYX)
    for i in 1:1000
        r = rot_type(100.0 * randn(), 100.0 * randn(), 100.0 * randn())
        r_prin = principal_value(r)
        @test -π ≤ r_prin.theta1 ≤ π
        @test -π ≤ r_prin.theta2 ≤ π
        @test -π ≤ r_prin.theta3 ≤ π
        @test r_prin ≈ r
    end
end


text:
import Base.show

export BinaryFileHeader, fh_byte2sample, show

mutable struct BinaryFileHeader
    Job                                  :: Int32
    Line                                 :: Int32
    Reel                                 :: Int32
    DataTracePerEnsemble                 :: Int16
    AuxiliaryTracePerEnsemble            :: Int16
    dt                                   :: Int16
    dtOrig                               :: Int16
    ns                                   :: Int16
    nsOrig                               :: Int16
    DataSampleFormat                     :: Int16
    EnsembleFold                         :: Int16
    TraceSorting                         :: Int16
    VerticalSumCode                      :: Int16
    SweepFrequencyStart                  :: Int16
    SweepFrequencyEnd                    :: Int16
    SweepLength                          :: Int16
    SweepType                            :: Int16
    SweepChannel                         :: Int16
    SweepTaperlengthStart                :: Int16
    SweepTaperLengthEnd                  :: Int16
    TaperType                            :: Int16
    CorrelatedDataTraces                 :: Int16
    BinaryGain                           :: Int16
    AmplitudeRecoveryMethod              :: Int16
    MeasurementSystem                    :: Int16
    ImpulseSignalPolarity                :: Int16
    VibratoryPolarityCode                :: Int16
    SegyFormatRevisionNumber             :: Int16
    FixedLengthTraceFlag                 :: Int16
    NumberOfExtTextualHeaders            :: Int16
end

function BinaryFileHeader()

   BinaryFileHeader(0,0,0,0,0,0,0,0,0,0,
                    0,0,0,0,0,0,0,0,0,0,
                    0,0,0,0,0,0,0,0,0,0)

end

function fh_byte2sample()

    Dict{String, Int32}(
        "Job"                               => 3200,
        "Line"                              => 3204,
        "Reel"                              => 3208,
        "DataTracePerEnsemble"              => 3212,
        "AuxiliaryTracePerEnsemble"         => 3214,
        "dt"                                => 3216,
        "dtOrig"                            => 3218,
        "ns"                                => 3220,
        "nsOrig"                            => 3222,
        "DataSampleFormat"                  => 3224,
        "EnsembleFold"                      => 3226,
        "TraceSorting"                      => 3228,
        "VerticalSumCode"                   => 3230,
        "SweepFrequencyStart"               => 3232,
        "SweepFrequencyEnd"                 => 3234,
        "SweepLength"                       => 3236,
        "SweepType"                         => 3238,
        "SweepChannel"                      => 3240,
        "SweepTaperlengthStart"             => 3242,
        "SweepTaperLengthEnd"               => 3244,
        "TaperType"                         => 3246,
        "CorrelatedDataTraces"              => 3248,
        "BinaryGain"                        => 3250,
        "AmplitudeRecoveryMethod"           => 3252,
        "MeasurementSystem"                 => 3254,
        "ImpulseSignalPolarity"             => 3256,
        "VibratoryPolarityCode"             => 3258,
        "SegyFormatRevisionNumber"          => 3500,
        "FixedLengthTraceFlag"              => 3502,
        "NumberOfExtTextualHeaders"         => 3504)
end

function show(io::IO, bfh::BinaryFileHeader)

    println("BinaryFileHeader:")

    for field in fieldnames(BinaryFileHeader)
        s = @sprintf "     %30s: %9d" String(field) getfield(bfh, field)
        println(s)
    end
    println("\n")
end


text:
using MsgPack
using Base.Test
using Compat

ck_pack(a, b) = pack(a) == b && unpack(b) == a

# nil
@test ck_pack(nothing, [0xc0])

# bool
@test ck_pack(false, [0xc2])
@test ck_pack(true, [0xc3])

# positive fixint
@test ck_pack(5, [0x05])
# UInt8
@test ck_pack(128, [0xcc,0x80])
# UInt16
@test ck_pack(2^8, [0xcd,0x01,0x00])
# UInt32
@test ck_pack(2^16, [0xce,0x00,0x01,0x00,0x00])
# UInt64
@test ck_pack(2^32, [0xcf,0x00,0x00,0x00,0x01,0x00,0x00,0x00,0x00])
@test ck_pack((@compat UInt64(2)^64-1), [0xcf,0xff,0xff,0xff,0xff,0xff,0xff,0xff,0xff])

# negative fixint
@test ck_pack(-5, [0xfb])
@test ck_pack(-32, [0xe0])
# Int8
@test ck_pack(-33, [0xd0,0xdf])
@test ck_pack(-2^7, [0xd0,0x80])
# Int16
@test ck_pack(-2^8, [0xd1,0xff,0x00])
@test ck_pack(-2^15, [0xd1,0x80,0x00])
# Int32
@test ck_pack(-2^16, [0xd2,0xff,0xff,0x00,0x00])
@test ck_pack(-2^31, [0xd2,0x80,0x00,0x00,0x00])
# Int64
@test ck_pack(-2^32, [0xd3,0xff,0xff,0xff,0xff,0x00,0x00,0x00,0x00])
@test ck_pack(-2^63, [0xd3,0x80,0x00,0x00,0x00,0x00,0x00,0x00,0x00])

# Float32
@test ck_pack((@compat Float32(5.876)), [0xca, 0x40, 0xbc, 0x08, 0x31])
# Float64
@test ck_pack(0.654789321, [0xcb, 0x3f, 0xe4, 0xf4, 0x08, 0xbb, 0xee, 0xe1, 0xa8])

# fixstr
@test ck_pack("Меня зовут Шон",
              [0xba,0xd0,0x9c,0xd0,0xb5,0xd0,0xbd,0xd1,0x8f,
               0x20,0xd0,0xb7,0xd0,0xbe,0xd0,0xb2,0xd1,0x83,
               0xd1,0x82,0x20,0xd0,0xa8,0xd0,0xbe,0xd0,0xbd])
# str8 - currently unimplemented
# str16
s = randstring(2^8)
@test ck_pack(s, vcat(0xda, 0x01, 0x00, convert(Vector{UInt8}, s)))
s = randstring(2^16 - 1)
@test ck_pack(s, vcat(0xda, 0xff, 0xff, convert(Vector{UInt8}, s)))
# str32
s = randstring(2^16)
@test ck_pack(s, vcat(0xdb, 0x00, 0x01, 0x00, 0x00, convert(Vector{UInt8}, s)))
s = ""

# bin8
b = rand(UInt8, 121)
@test ck_pack(b, vcat(0xc4, 0x79, b))
# bin16
b = rand(UInt8, 2^11)
@test ck_pack(b, vcat(0xc5, 0x08, 0x00, b))
# bin32
b = rand(UInt8, 2^20)
@test ck_pack(b, vcat(0xc6, 0x00, 0x10, 0x00, 0x00, b))
b = []

# fixarray
@test ck_pack(Any[1, (@compat Float32(2)), UInt8[0x22, 0x54]],
              [0x93, 0x01, 0xca, 0x40, 0x00, 0x00, 0x00, 0xc4, 0x02, 0x22, 0x54])
@test ck_pack(Any[nothing, "heya", true, false],
              [0x94, 0xc0, 0xa4, 0x68, 0x65, 0x79, 0x61, 0xc3, 0xc2])
@test ck_pack(Any[true, UInt8[0xff, 0xde, 0x11], [1, 2, 3]],
              [0x93, 0xc3, 0xc4, 0x03, 0xff, 0xde, 0x11, 0x93, 0x01, 0x02, 0x03])

# tuple
@test pack(Any[(1, [true, false], UInt8[0xff, 0xde])]) ==
              [0x91, 0x93, 0x01, 0x92, 0xc3, 0xc2, 0xc4, 0x02, 0xff, 0xde]
@test pack(("hi", ((@compat Float32(2)), 0xff))) ==
              [0x92, 0xa2, 0x68, 0x69, 0x92, 0xca, 0x40, 0x00, 0x00, 0x00, 0xcc, 0xff]

# fixmap
@test ck_pack((@compat Dict(1=>2, "hi"=>"mom")),
              [0x82,0xa2,0x68,0x69,0xa3,0x6d,0x6f,0x6d,0x01,0x02])

# fixext 1
@test ck_pack(Ext(4, [0xbb]), [0xd4, 0x04, 0xbb])
# fixext 2
@test ck_pack(Ext(-32, [0x56, 0x8d], impltype=true), [0xd5, 0xe0, 0x56, 0x8d])
# fixext 4
@test ck_pack(Ext(-123, [0x80, 0x7c, 0x8b, 0xf8], impltype=true),
              [0xd6, 0x85, 0x80, 0x7c, 0x8b, 0xf8])
# fixext 8
@test ck_pack(Ext(111, [0x04, 0x16, 0x94, 0x13, 0x0a, 0x7d, 0x6f, 0x0c]),
              [0xd7, 0x6f, 0x04, 0x16, 0x94, 0x13, 0x0a, 0x7d, 0x6f, 0x0c])
# fixext 16
@test ck_pack(Ext(79, [0x00, 0x30, 0xd5, 0x64, 0x0f, 0x8d, 0x92, 0x90,
                       0x98, 0x99, 0x14, 0x57, 0x0e, 0x8d, 0xf1, 0x3a]),
              [0xd8, 0x4f, 0x00, 0x30, 0xd5, 0x64, 0x0f, 0x8d, 0x92,
               0x90, 0x98, 0x99, 0x14, 0x57, 0x0e, 0x8d, 0xf1, 0x3a])
# no elements
@test ck_pack(Ext(-118, UInt8[], impltype=true), [0xc7, 0x00, 0x8a])
# ext 8
@test ck_pack(Ext(50, [0x62, 0x4c, 0x7c, 0x0f, 0x86, 0x04]),
              [0xc7, 0x06, 0x32, 0x62, 0x4c, 0x7c, 0x0f, 0x86, 0x04])
# ext 16
b = rand(UInt8, 2^14)
@test ck_pack(Ext(78, b), vcat(0xc8, 0x40, 0x00, 0x4e, b))
# ext 32
b = rand(UInt8, 2^19)
@test ck_pack(Ext(-123, b, impltype=true),
              vcat(0xc9, 0x00, 0x08, 0x00, 0x00, 0x85, b))


text:
using HiQGA.transD_GP, PyPlot, DelimitedFiles, DataInterpolations
## model fixed parts, i.e., air, but only to be compatible with AEM ...
zfixed   = [-1e5]
ρfixed   = [1e12]
nmax = 200
# z grid spec starts, the first z and first ρ will be unused in MT
zstart = 0.0
extendfrac, dz = 1.169, 10
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=50, showplot=true)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
## SERPENT S08
D = readdlm("MTS08.txt")
T = 1 ./D[:,1]
d_log10_ρ = D[:,2]
d_phase_deg = -D[:,3]
σ_log10_ρ = D[:,4]
σ_phase_deg = D[:,5]
## load zrho limits
usedepthprior = true # toggle this and see differences in posterior
zrholim = readdlm("zlog10rholims.txt")
interplow, interphigh = map(x->ConstantInterpolation(x, zrholim[:,1]), (zrholim[:,2], zrholim[:,3]))
ρlow, ρhigh = interplow.(zall), interphigh.(zall)
Δ = ρhigh - ρlow
F = transD_GP.MT1DInversion.MT1DZ_nodepthprior(d_log10_ρ = d_log10_ρ, d_phase_deg = d_phase_deg,
                                    σ_log10_ρ = σ_log10_ρ, σ_phase_deg = σ_phase_deg,
                                    freqs = 1 ./T, zboundaries = zboundaries, irxlayer = 1)
F = transD_GP.MT1DInversion.makestretchop(F, ρlow=ρlow, Δ=Δ) # apply the stretch prior and get new operator
transD_GP.MT1DInversion.plotdata(F)
F.stretch = usedepthprior
if usedepthprior
    transD_GP.MT1DInversion.plotpriorenv(F, lc = "r", plotlinear=false)
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
# include(Pkg.dir("KUnet/test/mnist.jl"))

reload("src/NNN.jl")

const batch_size  = 784
const input_size  = 1
const output_size = 10
const hidden_size = 100

function toy_data(n)
    x = zeros(Float32, input_size,  n*batch_size)
    y = zeros(Float32, output_size, n*batch_size)

    for i in 0:batch_size:((n-1)*batch_size)
        mnis_i = floor(Int64, rand() * 50000) + 1
        mnis_x = MNIST.xtrn[:, mnis_i]
        mnis_y = MNIST.ytrn[:, mnis_i]

        for j in 1:batch_size
            pos = i + j
            x[:, pos] = mnis_x[j]
            y[:, pos] = mnis_y
        end
    end

    (x, y)
end

function rnn_train(xtrn, ytrn, xtst, ytst, iter)
    net = NNN.Network(batch_size)
    NNN.push!(net,
             NNN.SimpleRecurrentLayer(batch_size, input_size, hidden_size;
                                     activation=NNN.LeakyRelu()),
             NNN.AdaGrad(input_size, hidden_size;
                        ε=0.2,
                        μ=0.05,
                        λ=0.0000001),
             NNN.AdaGrad(hidden_size, hidden_size;
                        ε=0.2,
                        μ=0.05,
                        λ=0.000006))
    NNN.push!(net,
             NNN.OutputLayer(batch_size, hidden_size, output_size;
                            activation=NNN.Softmax(),
                            loss=NNN.rnn_cross_entropy_cost),
             NNN.SGD(;λ=0.000001,
                    ε=0.05,
                    gradient_clip=10.0f0))

    @time NNN.train(xtrn, ytrn, net; iter=iter, batch_size=batch_size)
    errors = NNN.test(xtst, ytst, net; batch_size=batch_size)

    (net, errors)
end

xtrn, ytrn = toy_data(2000);
xtst, ytst = toy_data(100);
net, errors = rnn_train(xtrn, ytrn, xtst, ytst, 1); errors

# NNN.train(xtrn, ytrn, net; iter=10, batch_size=batch_size)


text:
using SeisPlot, SeisAcoustic

# ==============================================================================
#             read physical model
dir_work = joinpath(homedir(), "Desktop");

path_rho = joinpath(dir_work, "physical_model/rho.rsf");
path_vel = joinpath(dir_work, "physical_model/vel.rsf");

# read the physical model
println("Loading velocity and density model");
(hdr_rho, rho) = read_RSdata(path_rho);
(hdr_vel, vel) = read_RSdata(path_vel);

# # cropped model for imaging
z_range = 1:3:300;
x_range = 4000:2:4600;
vel = vel[z_range, x_range];
rho = rho[z_range, x_range];
vel = model_smooth(vel, 10);

# used for removing direct arrival
rho1 = copy(rho); rho1 .= minimum(rho);
vel1 = copy(vel); vel1 .= vel[1];

# SeisPlotTX(vel, hbox=1.5*2, wbox=3*2, cmap="rainbow", vmin=minimum(vel), vmax=maximum(vel));
# SeisPlotTX(rho, hbox=1.5*2, wbox=3*2, cmap="rainbow", vmin=minimum(rho), vmax=maximum(rho));

# vertical and horizontal grid size
dz = 6.25; dx = 6.25;

# time step size and maximum modelling length
dt = 0.0008; tmax = 1.25;

# top absorbing boundary to aviod surface-related multiple
free_surface = false;
data_format  = Float32;
order        = 2;

# tdparams for generating observations
fidiff_hete = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
                       data_format=data_format, order=order);

# tdparams for removing direct arrival
fidiff_homo = TdParams(rho1, vel1, free_surface, dz, dx, dt, tmax;
                       data_format=data_format, order=order);

# tdparams for imaging
fidiff = TdParams(rho1, vel, free_surface, dz, dx, dt, tmax;
                  data_format=data_format, order=order);

# vector of source
println("Set acquisition geometry");
isx = collect(30 : 70 : fidiff.nx-30); isz = 2*ones(length(isx));
src = get_multi_sources(isz, isx, fidiff; amp=100000, fdom=20, type_flag="miniphase");

# receiver location
irx = collect(1: 2 : fidiff.nx);
irz = 2 * ones(Int64, length(irx));

dir_obs        = joinpath(dir_work, "observations");
dir_sourceside = joinpath(dir_work, "sourceside");

# simulate observed data (direct arrival removed)
println("Simulate observations");
get_reflections(dir_obs, irz, irx, src, fidiff_hete, fidiff_homo);

println("Compute source-side wavefield and save its boundary");
get_wavefield_bound(dir_sourceside, src, fidiff);

# dobs = read_recordings(joinpath(dir_obs, "recordings_1.bin"));
# SeisPlotTX(dobs.p, dy=dt, cmap="gray", wbox=6, hbox=6);

# lsrtm
println("CGLS");
born_params = (irz=irz, irx=irx, dir_sourceside=dir_sourceside, fidiff=fidiff, normalization_flag=true, mute_index=10);
(x, his) = cgls(born_approximation, dir_obs; dir_work=dir_work, op_params=born_params, maxIter=25,
                d_axpby=recordings_axpby!, m_axpby=image_axpby!, d_norm=recordings_norm, m_norm=l2norm_rsf);

# check the result
println("plotting");
path = joinpath(dir_work, "sourceside/normalization.rsf");
(hdr, s) = read_RSdata(path);

path = joinpath(dir_work, "iterations/iteration_1.rsf");
(hdr, m) = read_RSdata(path);
p = m .* vec(s);
p = reshape(p, fidiff.nz, fidiff.nx);
p1= laplace_filter(p);
SeisPlotTX(p1, cmap="gray", wbox=10, hbox=5, title="RTM", name=joinpath(dir_work, "rtm.pdf"));

path = joinpath(dir_work, "iterations/iteration_10.rsf");
(hdr, m) = read_RSdata(path);
p = m .* vec(s);
p = reshape(p, fidiff.nz, fidiff.nx);
p1= laplace_filter(p);
SeisPlotTX(p1, cmap="gray", wbox=10, hbox=5, title="LSRTM_10", name=joinpath(dir_work, "lsrtm_10.pdf"));

path = joinpath(dir_work, "iterations/iteration_25.rsf");
(hdr, m) = read_RSdata(path);
p = m .* vec(s);
p = reshape(p, fidiff.nz, fidiff.nx);
p1= laplace_filter(p);
SeisPlotTX(p1, cmap="gray", wbox=10, hbox=5, title="LSRTM_25", name=joinpath(dir_work, "lsrtm_25.pdf"));

plot(his); savefig(joinpath(dir_work, "convergence.pdf")); close();


text:
# workspace()
using Indicators
using PyPlot
using Dates
using Temporal

aapl = yahoo("AAPL")
aapl = aapl["2015"]

t = aapl.index
aapl = aapl.values
Open = aapl[:,1]
High = aapl[:,2]
Low = aapl[:,3]
Close = aapl[:,4]
Volume = aapl[:,5]
HLC = [High Low Close]

subplot(411)
plot(t, Close, lw=2, c="k", label="AAPL")
plot(t, kama(Close), c="b", label="Kaufman AMA")
plot(t, trima(Close), c="g", label="Triangula MA")
plot(t, hma(Close), c="r", label="Hull MA")
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

subplot(412)
plot(t, kst(Close), c="m", label="KST")
plot(t, sma(kst(Close),n=9), c="c", label="Signal")
plot([t[1],t[end]], [0,0], ls="--", c=[0.4,0.4,0.4])
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

subplot(413)
plot(t, wpr(HLC), c=[1,0.5,0], label="Williams %R")
plot([t[1],t[end]], [-20,-20], c="r", ls="--")
plot([t[1],t[end]], [-80,-80], c="g", ls="--")
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

subplot(414)
plot(t, cci(HLC), c="c", label="CCI")
plot([t[1],t[end]], [-100,-100], c="g", ls="--")
plot([t[1],t[end]], [100,100], c="r", ls="--")
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

tight_layout()


text:
using FillArrays, FFTW
using FillArrays: AbstractFill, getindex_value
using Base.Broadcast: broadcasted, broadcast_shape

@adjoint (::Type{T})(::UndefInitializer, args...) where T<:Array = T(undef, args...), Δ -> nothing

@adjoint Array(xs::AbstractArray) = Array(xs), ȳ -> (ȳ,)
@adjoint Array(xs::Array) = Array(xs), ȳ -> (ȳ,)

@nograd size, length, eachindex, axes, Colon(), findfirst, findlast, findall, randn, ones, zeros, one, zero,
  print, println, any, all

@adjoint rand(dims::Integer...) = rand(dims...), _ -> nothing

@adjoint Base.vect(xs...) = Base.vect(xs...), Δ -> (Δ...,)

@adjoint copy(x::AbstractArray) = copy(x), ȳ -> (ȳ,)

# Array Constructors
@adjoint (::Type{T})(x::T) where T<:Array = T(x), ȳ -> (ȳ,)
@adjoint function (::Type{T})(x::Number, sz) where {T <: Fill}
    back(Δ::AbstractArray) = (sum(Δ), nothing)
    back(Δ::NamedTuple) = (Δ.value, nothing)
    return Fill(x, sz), back
end

@adjoint (::Type{T})(sz) where {T<:Zeros} = Zeros(sz), Δ->(nothing,)
@adjoint (::Type{T})(sz) where {T<:Ones} = Ones(sz), Δ->(nothing,)

_zero(xs::AbstractArray{<:Number}, T=float(eltype(xs))) = fill!(similar(xs, T), false)
_zero(xs::AbstractArray, T=Any) = Union{Nothing, T}[nothing for x in xs]

@adjoint getindex(x::AbstractArray, inds...) = x[inds...], ∇getindex(x, inds)

@adjoint view(x::AbstractArray, inds...) = view(x, inds...), ∇getindex(x, inds)

∇getindex(x::AbstractArray, inds) = dy -> begin
  if inds isa  NTuple{<:Any,Integer}
    dx = _zero(x, typeof(dy))
    dx[inds...] = dy
  else
    dx = _zero(x, eltype(dy))
    @views dx[inds...] .+= dy
  end
  (dx, map(_->nothing, inds)...)
end

@adjoint! setindex!(xs::AbstractArray, x...) = setindex!(xs, x...),
  _ -> error("Mutating arrays is not supported")

for f in [push!, pop!, pushfirst!, popfirst!]
  @eval @adjoint! $f(xs::Vector, x...) =
    push!(xs, x...), _ -> error("Mutating arrays is not supported")
end

# General

@adjoint collect(x::Array) = collect(x), Δ -> (Δ,)

@adjoint fill(x::Real, dims...) = fill(x, dims...), Δ->(sum(Δ), map(_->nothing, dims)...)

@adjoint function circshift(A, shifts)
  circshift(A, shifts), Δ -> (circshift(Δ, map(-, shifts)), nothing)
end

@adjoint permutedims(xs) = permutedims(xs), Δ -> (permutedims(Δ),)

@adjoint permutedims(xs::AbstractVector) = permutedims(xs), Δ -> (vec(permutedims(Δ)),)

@adjoint permutedims(xs, dims) = permutedims(xs, dims),
  Δ -> (permutedims(Δ, invperm(dims)), nothing)

@adjoint PermutedDimsArray(xs, dims) = PermutedDimsArray(xs, dims),
  Δ -> (PermutedDimsArray(Δ, invperm(dims)), nothing)

@adjoint reshape(xs, dims...) = reshape(xs, dims...),
  Δ -> (reshape(Δ, size(xs)),map(_->nothing,dims)...)

@adjoint function hvcat(rows::Tuple{Vararg{Int}}, xs::T...) where T<:Number
  hvcat(rows, xs...), ȳ -> (nothing, permutedims(ȳ)...)
end

pull_block_vert(sz, Δ, A::AbstractVector) = Δ[sz-length(A)+1:sz]
pull_block_vert(sz, Δ, A::AbstractMatrix) = Δ[sz-size(A, 1)+1:sz, :]
@adjoint function vcat(A::Union{AbstractVector, AbstractMatrix}...)
  sz = cumsum([size.(A, 1)...])
  return vcat(A...), Δ->(map(n->pull_block_vert(sz[n], Δ, A[n]), eachindex(A))...,)
end
@adjoint vcat(xs::Number...) = vcat(xs...), Δ -> (Δ...,)

pull_block_horz(sz, Δ, A::AbstractVector) = Δ[:, sz]
pull_block_horz(sz, Δ, A::AbstractMatrix) = Δ[:, sz-size(A, 2)+1:sz]
@adjoint function hcat(A::Union{AbstractVector, AbstractMatrix}...)
  sz = cumsum([size.(A, 2)...])
  return hcat(A...), Δ->(map(n->pull_block_horz(sz[n], Δ, A[n]), eachindex(A))...,)
end
@adjoint hcat(xs::Number...) = hcat(xs...), Δ -> (Δ...,)

@adjoint function cat(Xs...; dims)
  cat(Xs...; dims = dims), Δ -> begin
    start = ntuple(_ -> 0, ndims(Δ))
    dXs = map(Xs) do x
      move = ntuple(d -> d in dims ? size(x,d) : 0, ndims(Δ))
      x_in_Δ = ntuple(d -> move[d] > 0 ? (start[d]+1:start[d]+move[d]) : Colon(), ndims(Δ))
      start = start .+ move
      dx = reshape(Δ[x_in_Δ...], size(x))
    end
  end
end

@adjoint function repeat(xs; inner=ntuple(_->1, ndims(xs)), outer=ntuple(_->1, ndims(xs)))
  repeat(xs, inner = inner, outer = outer), function (Δ)
    Δ′ = zero(xs)
    S = size(xs)

    # Loop through each element of Δ, calculate source dimensions, accumulate into Δ′
    for (dest_idx, val) in pairs(IndexCartesian(), Δ)
        # First, round dest_idx[dim] to nearest gridpoint defined by inner[dim], then
        # wrap around based on original size S.
        src_idx = [mod1(div(dest_idx[dim] - 1, inner[dim]) + 1, S[dim]) for dim in 1:length(S)]
        Δ′[src_idx...] += val
    end
    return (Δ′,)
  end
end

@adjoint getindex(i::Int, j::Int) = i[j], _ -> nothing

function unzip(tuples)
  map(1:length(first(tuples))) do i
      map(tuple -> tuple[i], tuples)
  end
end
function ∇map(cx, f, args...)
  ys_and_backs = map((args...) -> _pullback(cx, f, args...), args...)
  if isempty(ys_and_backs)
    ys_and_backs, _ -> nothing
  else
    ys, backs = unzip(ys_and_backs)
    ys, function (Δ)
      Δf_and_args_zipped = map((f, δ) -> f(δ), backs, Δ)
      Δf_and_args = unzip(Δf_and_args_zipped)
      Δf = reduce(accum, Δf_and_args[1])
      (Δf, Δf_and_args[2:end]...)
    end
  end
end

@adjoint function map(f, args::Union{AbstractArray,Tuple}...)
  ∇map(__context__, f, args...)
end

function _pullback(cx::AContext, ::typeof(collect), g::Base.Generator)
  y, back = ∇map(cx, g.f, g.iter)
  y, function (ȳ)
    f̄, x̄ = back(ȳ)
    (nothing, (f = f̄, iter = x̄),)
  end
end

@adjoint iterate(r::UnitRange, i...) = iterate(r, i...), _ -> nothing

@adjoint function sort(x::AbstractArray)
  p = sortperm(x)
  return x[p], x̄ -> (x̄[invperm(p)],)
end

# Reductions

@adjoint function sum(xs::AbstractArray; dims = :)
  if dims === (:)
    sum(xs), Δ -> (Fill(Δ, size(xs)),)
  else
    sum(xs, dims = dims), Δ -> (similar(xs) .= Δ,)
  end
end

function _pullback(cx::AContext, ::typeof(sum), f, xs::AbstractArray)
  y, back = pullback(cx, (xs -> sum(f.(xs))), xs)
  y, ȳ -> (nothing, nothing, back(ȳ)...)
end

@adjoint function sum(::typeof(abs2), X::AbstractArray; dims = :)
  return sum(abs2, X; dims=dims), Δ::Union{Number, AbstractArray}->(nothing, ((2Δ) .* X))
end

@adjoint function prod(xs::AbstractArray{<:Number}; dims = :)
  if dims === (:)
    prod(xs), Δ -> (prod(xs) ./ xs .* Δ,)
  else
    prod(xs, dims = dims),
      Δ -> (reshape(.*(circshift.([reshape(xs, length(xs))], 1:length(xs)-1)...), size(xs)) .* Δ,)
  end
end

function _pullback(cx::AContext, ::typeof(prod), f, xs::AbstractArray)
  y, back = pullback(cx, (xs -> prod(f.(xs))), xs)
  y, ȳ -> (nothing, nothing, back(ȳ)...)
end

@adjoint function maximum(xs::AbstractArray; dims = :)
  max, i = findmax(xs, dims = dims)
  max, function (Δ)
    Δ isa Real && abs(Δ) <= sqrt(eps(float(Δ))) && return nothing
    Δ′ = zero(xs)
    Δ′[i] = Δ
    return (Δ′,)
  end
end

@adjoint function minimum(xs::AbstractArray; dims = :)
  min, i = findmin(xs, dims = dims)
  min, function (Δ)
    Δ′ = zero(xs)
    Δ′[i] = Δ
    return (Δ′,)
  end
end

@adjoint function dropdims(xs::AbstractArray; dims)
  dropdims(xs, dims = dims), Δ -> (reshape(Δ, size(xs)...),)
end

@adjoint real(x::AbstractArray) = real(x), r̄ -> (real(r̄),)
@adjoint conj(x::AbstractArray) = conj(x), r̄ -> (conj(r̄),)
@adjoint imag(x::AbstractArray) = imag(x), ī -> (complex.(0, real.(ī)),)

@adjoint function mean(xs::AbstractArray; dims = :)
  return mean(xs, dims=dims), Δ -> (_backmean(xs,Δ,dims),)
end
_backmean(xs, Δ, ::Colon) = zero(xs) .+ Δ ./ length(xs)
_backmean(xs, Δ, dims) = zero(xs) .+ Δ ./ mapreduce(i -> size(xs,i),*,dims)

@adjoint function Statistics.var(xs::AbstractArray; corrected::Bool=true, dims=:, mean=mean(xs, dims=dims))
    return Statistics.var(xs; corrected=corrected, mean=mean, dims=dims), Δ -> _backvar(xs, Δ, corrected, mean, dims)
end
_backvar(xs, Δ, corrected::Bool, mean, dims)         = _backvar(xs, Δ, mapreduce(i -> size(xs,i),*,dims) - corrected, mean)
_backvar(xs, Δ, corrected::Bool, mean, ::Colon)      = _backvar(xs, Δ, length(xs) - corrected, mean)
_backvar(xs, Δ, N::Int, mean) = (convert(eltype(xs), 2/N) .* Δ .* (xs .- mean),)

@adjoint function Statistics.std(xs::AbstractArray; corrected::Bool=true, dims=:, mean=mean(xs, dims=dims))
    s = Statistics.std(xs; corrected=corrected, mean=mean, dims=dims)
    return s, Δ -> _backvar(xs, Δ ./ (2 .* s), corrected, mean, dims)
end

# LinAlg
# ======

@adjoint function(A::AbstractMatrix * B::AbstractMatrix)
  return A * B, Δ::AbstractMatrix->(Δ * B', A' * Δ)
end

@adjoint function(A::AbstractMatrix * x::AbstractVector)
  return A * x, Δ::AbstractVector->(Δ * x', A' * Δ)
end

@adjoint function *(x::Union{Transpose{<:Any, <:AbstractVector},
                             LinearAlgebra.Adjoint{<:Any, <:AbstractVector}},
                    y::AbstractVector)
  return x * y, Δ->(Δ * y', x' * Δ)
end

@adjoint function(a::AbstractVector * x::AbstractMatrix)
  return a * x, Δ::AbstractMatrix->(vec(Δ * x'), a' * Δ)
end

@adjoint function transpose(x)
  back(Δ) = (transpose(Δ),)
  back(Δ::NamedTuple{(:parent,)}) = (Δ.parent,)
  return transpose(x), back
end

@adjoint function Base.adjoint(x)
  back(Δ) = (Δ',)
  back(Δ::NamedTuple{(:parent,)}) = (Δ.parent,)
  return x', back
end

@adjoint parent(x::LinearAlgebra.Adjoint) = parent(x), ȳ -> (LinearAlgebra.Adjoint(ȳ),)

@adjoint dot(x::AbstractArray, y::AbstractArray) = dot(x, y), Δ->(Δ .* y, Δ .* x)

function _kron(mat1::AbstractMatrix,mat2::AbstractMatrix)
    m1, n1 = size(mat1)
    mat1_rsh = reshape(mat1,(1,m1,1,n1))

    m2, n2 = size(mat2)
    mat2_rsh = reshape(mat2,(m2,1,n2,1))

    return reshape(mat1_rsh.*mat2_rsh, (m1*m2,n1*n2))
end

@adjoint kron(a::AbstractMatrix, b::AbstractMatrix) = pullback(_kron, a, b)

@adjoint function Diagonal(d::AbstractVector)
  back(Δ::NamedTuple) = (Δ.diag,)
  back(Δ::AbstractMatrix) = (diag(Δ),)
  return Diagonal(d), back
end

@adjoint diag(A::AbstractMatrix) = diag(A), Δ->(Diagonal(Δ),)

@adjoint det(xs) = det(xs), Δ -> (Δ * det(xs) * transpose(inv(xs)),)

@adjoint logdet(xs) = logdet(xs), Δ -> (Δ * transpose(inv(xs)),)

@adjoint logabsdet(xs) = logabsdet(xs), Δ -> (Δ[1] * transpose(inv(xs)),)

@adjoint function inv(A)
  Ainv = inv(A)
  return Ainv, function (Δ)
    ∇A = - Ainv' * Δ * Ainv'
    return (∇A, )
  end
end

# Defaults for atol and rtol copied directly from LinearAlgebra. See the following for
# derivation:
# Golub, Gene H., and Victor Pereyra. "The differentiation of pseudo-inverses and nonlinear
# least squares problems whose variables separate." SIAM Journal on numerical analysis 10.2
# (1973): 413-432.
@adjoint function pinv(
  A::AbstractMatrix{T};
  atol::Real = 0.0,
  rtol::Real = (eps(real(float(one(T))))*min(size(A)...))*iszero(atol),
) where {T}
  Y = pinv(A)
  return Y, Δ->(-Y' * Δ * Y' + (I - A * Y) * Δ' * Y * Y' + Y' * Y * Δ' * (I - Y * A),)
end

# When `A` is guaranteed to be square, definitely use the simple expression for the adjoint.
@adjoint function \(
  A::Union{
    Diagonal,
    AbstractTriangular,
    LinearAlgebra.Adjoint{<:Any, <:AbstractTriangular},
    Transpose{<:Any, <:AbstractTriangular},
  },
  B::AbstractVecOrMat,
)
  Y = A \ B
  return Y, function(Ȳ)
    B̄ = A' \ Ȳ
    return (-B̄ * Y', B̄)
  end
end

@adjoint function /(A::AbstractMatrix, B::Union{Diagonal, AbstractTriangular})
  Y = A / B
  return Y, function(Ȳ)
    Ā = Ȳ / B'
    return (Ā, -Y' * Ā)
  end
end

@adjoint function \(A::AbstractMatrix, B::AbstractVecOrMat)
  Z = A \ B
  return Z, function(Z̄)
    B̄ = A' \ Z̄
    if size(A, 1) == size(A, 2)
      return (-B̄ * Z', B̄)
    else
      a = -B̄ * Z'
      b = (B - A * Z) * B̄' / A'
      c = A' \ Z * (Z̄' - B̄' * A)
      return (a + b + c, B̄)
    end
  end
end

function _pullback(cx::AContext, ::typeof(norm), x::AbstractArray, p::Real = 2)
  fallback = (x, p) -> sum(abs.(x).^p .+ eps(0f0))^(1/p) # avoid d(sqrt(x))/dx == Inf at 0
  _pullback(cx, fallback, x, p)
end

# LinAlg Matrix Types
# ===================

@adjoint LinearAlgebra.LowerTriangular(A) = LowerTriangular(A), Δ->(LowerTriangular(Δ),)
@adjoint LinearAlgebra.UpperTriangular(A) = UpperTriangular(A), Δ->(UpperTriangular(Δ),)
@adjoint LinearAlgebra.UnitLowerTriangular(A) = UnitLowerTriangular(A), Δ->(UnitLowerTriangular(Δ)-I,)
@adjoint LinearAlgebra.UnitUpperTriangular(A) = UnitUpperTriangular(A), Δ->(UnitUpperTriangular(Δ)-I,)

# This is basically a hack while we don't have a working `ldiv!`.
@adjoint function \(A::Cholesky, B::AbstractVecOrMat)
  Y, back = Zygote.pullback((U, B)->U \ (U' \ B), A.U, B)
  return Y, function(Ȳ)
    Ā_factors, B̄ = back(Ȳ)
    return ((uplo=nothing, status=nothing, factors=Ā_factors), B̄)
  end
end

function _symmetric_back(Δ, uplo)
  L, U, D = LowerTriangular(Δ), UpperTriangular(Δ), Diagonal(Δ)
  return uplo == 'U' ? U .+ transpose(L) - D : L .+ transpose(U) - D
end
_symmetric_back(Δ::Diagonal, uplo) = Δ
_symmetric_back(Δ::UpperTriangular, uplo) = collect(uplo == 'U' ? Δ : transpose(Δ))
_symmetric_back(Δ::LowerTriangular, uplo) = collect(uplo == 'U' ? transpose(Δ) : Δ)

@adjoint function Symmetric(A::AbstractMatrix, uplo=:U)
  S = Symmetric(A, uplo)
  back(Δ::AbstractMatrix) = (_symmetric_back(Δ, S.uplo), nothing)
  back(Δ::NamedTuple) = (_symmetric_back(Δ.data, S.uplo), nothing)
  return S, back
end

_extract_imag(x) = (x->complex(0, imag(x))).(x)
function _hermitian_back(Δ, uplo)
  isreal(Δ) && return _symmetric_back(Δ, uplo)
  L, U, rD = LowerTriangular(Δ), UpperTriangular(Δ), real.(Diagonal(Δ))
  return uplo == 'U' ? U .+ L' - rD : L .+ U' - rD
end
_hermitian_back(Δ::Diagonal, uplo) = real.(Δ)
function _hermitian_back(Δ::LinearAlgebra.AbstractTriangular, uplo)
  isreal(Δ) && return _symmetric_back(Δ, uplo)
  ŪL̄ = Δ .- Diagonal(_extract_imag(diag(Δ)))
  if istriu(Δ)
    return collect(uplo == 'U' ? ŪL̄ : ŪL̄')
  else
    return collect(uplo == 'U' ? ŪL̄' : ŪL̄)
  end
end

@adjoint function LinearAlgebra.Hermitian(A::AbstractMatrix, uplo=:U)
  H = Hermitian(A, uplo)
  back(Δ::AbstractMatrix) = (_hermitian_back(Δ, H.uplo), nothing)
  back(Δ::NamedTuple) = (_hermitian_back(Δ.data, H.uplo), nothing)
  return H, back
end

@adjoint convert(::Type{R}, A::LinearAlgebra.HermOrSym{T,S}) where {T,S,R<:Array} = convert(R, A),
  Δ -> (nothing, convert(S, Δ),)
@adjoint Matrix(A::LinearAlgebra.HermOrSym{T,S}) where {T,S} = Matrix(A),
  Δ -> (convert(S, Δ),)

@adjoint function cholesky(Σ::Real)
  C = cholesky(Σ)
  return C, Δ::NamedTuple->(Δ.factors[1, 1] / (2 * C.U[1, 1]),)
end

@adjoint function cholesky(Σ::Diagonal)
  C = cholesky(Σ)
  return C, Δ::NamedTuple->(Diagonal(diag(Δ.factors) .* inv.(2 .* C.factors.diag)),)
end

# Implementation due to Seeger, Matthias, et al. "Auto-differentiating linear algebra."
@adjoint function cholesky(Σ::Union{StridedMatrix, Symmetric{<:Real, <:StridedMatrix}})
  C = cholesky(Σ)
  return C, function(Δ::NamedTuple)
    U, Ū = C.U, Δ.factors
    Σ̄ = Ū * U'
    Σ̄ = copytri!(Σ̄, 'U')
    Σ̄ = ldiv!(U, Σ̄)
    BLAS.trsm!('R', 'U', 'T', 'N', one(eltype(Σ)), U.data, Σ̄)
    @inbounds for n in diagind(Σ̄)
      Σ̄[n] /= 2
    end
    return (UpperTriangular(Σ̄),)
  end
end

@adjoint function lyap(A::AbstractMatrix, C::AbstractMatrix)
  X = lyap(A, C)
  return X, function (X̄)
    C̄ = lyap(collect(A'), X̄)
    Ā = C̄*X' + C̄'*X
    return (Ā, C̄)
  end
end

# Matrix of pairwise difference quotients
Base.@propagate_inbounds function _pairdiffquot(f, i, j, x, fx, dfx, d²fx = nothing)
  i == j && return dfx[i]
  Δx = x[i] - x[j]
  T = real(eltype(x))
  if d²fx === nothing
    abs(Δx) ≤ sqrt(eps(T)) && return (dfx[i] + dfx[j]) / 2
  else
    abs(Δx) ≤ eps(T)^(1/3) && return dfx[i] - Δx / 2 * d²fx[i]
  end
  Δfx = fx[i] - fx[j]
  return Δfx / Δx
end

Base.@propagate_inbounds function _pairdiffquotmat(f, n, x, fx, dfx, d²fx = nothing)
  Δfij = (i, j)->_pairdiffquot(f, i, j, x, fx, dfx, d²fx)
  return Δfij.(Base.OneTo(n), Base.OneTo(n)')
end

# Adjoint based on the Theano implementation, which uses the differential as described
# in Brančík, "Matlab programs for matrix exponential function derivative evaluation"
@adjoint exp(A::AbstractMatrix) = exp(A), function(F̄)
  n = size(A, 1)
  E = eigen(A)
  w = E.values
  ew = exp.(w)
  X = _pairdiffquotmat(exp, n, w, ew, ew, ew)
  V = E.vectors
  VF = factorize(V)
  Āc = (V * ((VF \ F̄' * V) .* X) / VF)'
  Ā = isreal(A) && isreal(F̄) ? real(Āc) : Āc
  return (Ā,)
end

@adjoint function LinearAlgebra.eigen(A::LinearAlgebra.RealHermSymComplexHerm)
  dU = eigen(A)
  return dU, function (Δ)
    d, U = dU
    d̄, Ū = Δ
    if Ū === nothing
      P = Diagonal(d̄)
    else
      F = inv.(d' .- d)
      P = F .* (U' * Ū)
      if d̄ === nothing
        P[diagind(P)] .= 0
      else
        P[diagind(P)] = d̄
      end
    end
    return (U * P * U',)
  end
end

@adjoint function LinearAlgebra.eigvals(A::LinearAlgebra.RealHermSymComplexHerm)
  d, U = eigen(A)
  return d, d̄ -> (U * Diagonal(d̄) * U',)
end

# Hermitian/Symmetric matrix functions that can be written as power series
_realifydiag!(A::AbstractArray{<:Real}) = A
function _realifydiag!(A)
  n = LinearAlgebra.checksquare(A)
  for i in 1:n
      @inbounds A[i,i] = real(A[i,i])
  end
  return A
end
@adjoint _realifydiag!(A) = _realifydiag!(A), Δ -> (_realifydiag!(Δ),)

_hasrealdomain(f, x) = true
_hasrealdomain(::Union{typeof.((acos,asin))...}, x) = all(x -> -1 ≤ x ≤ 1, x)
_hasrealdomain(::typeof(acosh), x) = all(x -> x ≥ 1, x)
_hasrealdomain(::Union{typeof.((log,sqrt,^))...}, x) = all(x -> x ≥ 0, x)

_process_series_eigvals(f, λ) = _hasrealdomain(f, λ) ? λ : complex.(λ)

_process_series_matrix(f, fA, A, fλ) = fA
_process_series_matrix(f, fA, ::LinearAlgebra.HermOrSym{<:Real}, fλ) = Symmetric(fA)
_process_series_matrix(f, fA, ::Hermitian{<:Complex}, ::AbstractVector{<:Real}) =
  Hermitian(_realifydiag!(fA))
_process_series_matrix(::typeof(^), fA, ::Hermitian{<:Real}, fλ) = Hermitian(fA)
_process_series_matrix(::typeof(^), fA, ::Hermitian{<:Real}, ::AbstractVector{<:Complex}) = fA
_process_series_matrix(::typeof(^), fA, ::Hermitian{<:Complex}, ::AbstractVector{<:Complex}) = fA

# Compute function on eigvals, thunks for conjugates of 1st and 2nd derivatives,
# and function to pull back adjoints to args
function _pullback_series_func_scalar(f, λ, args...)
  compλ = _process_series_eigvals(f, λ)
  fλ, fback = Zygote.pullback((x,args...) -> f.(x, args...), compλ, args...)
  n = length(λ)
  return (fλ,
          ()->fback(ones(n))[1],
          ()->nothing, # TODO: add 2nd deriv
          isempty(args) ? _ -> () : f̄λ -> tail(fback(f̄λ)))
end

function _pullback_series_func_scalar(f::typeof(^), λ, p)
  compλ = _process_series_eigvals(f, λ)
  r, powλ = isinteger(p) ? (Integer(p), λ) : (p, compλ)
  fλ = powλ .^ r
  return (fλ,
          ()->conj.(r .* powλ .^ (r - 1)),
          ()->conj.((r * (r - 1)) .* powλ .^ (r - 2)),
          f̄λ -> (dot(fλ .* log.(compλ), f̄λ),))
end

function _pullback_series_func_scalar(f::typeof(exp), λ)
  expλ = exp.(λ)
  return expλ, ()->expλ, ()->expλ, _ -> ()
end

_apply_series_func(f, A, args...) = f(A, args...)

@adjoint function _apply_series_func(f, A, args...)
  hasargs = !isempty(args)
  n = LinearAlgebra.checksquare(A)
  λ, U = eigen(A)
  fλ, dfthunk, d²fthunk, argsback = _pullback_series_func_scalar(f, λ, args...)
  fΛ = Diagonal(fλ)
  fA = U * fΛ * U'
  Ω = _process_series_matrix(f, fA, A, fλ)
  return Ω, function (f̄A)
    f̄Λ = U' * f̄A * U
    ārgs = hasargs ? argsback(diag(f̄Λ)) : ()
    P = _pairdiffquotmat(f, n, λ, conj(fλ), dfthunk(), d²fthunk())
    Ā = U * (P .* f̄Λ) * U'
    return (nothing, Ā, ārgs...)
  end
end

_hermsympow(A::Symmetric, p::Integer) = LinearAlgebra.sympow(A, p)
_hermsympow(A::Hermitian, p::Integer) = A^p

@adjoint function _hermsympow(A::Hermitian, p::Integer)
  if p < 0
    B, back = Zygote.pullback(A->Base.power_by_squaring(inv(A), -p), A)
  else
    B, back = Zygote.pullback(A->Base.power_by_squaring(A, p), A)
  end
  Ω = Hermitian(_realifydiag!(B))
  return Ω, function (Ω̄)
    B̄ = _hermitian_back(Ω̄, 'U')
    Ā = back(B̄)[1]
    return (Ā, nothing)
  end
end

_pullback(cx::AContext, ::typeof(^), A::LinearAlgebra.HermOrSym{<:Real}, p::Integer) =
  _pullback(cx, _hermsympow, A, p)
_pullback(cx::AContext, ::typeof(^), A::Symmetric{<:Complex}, p::Integer) =
  _pullback(cx, _hermsympow, A, p)
_pullback(cx::AContext, ::typeof(^), A::Hermitian{<:Complex}, p::Integer) =
  _pullback(cx, _hermsympow, A, p)

function _pullback(cx::AContext,
                   f::typeof(^),
                   A::LinearAlgebra.RealHermSymComplexHerm,
                   p::Real)
  return _pullback(cx, (A, p) -> _apply_series_func(f, A, p), A, p)
end

for func in (:exp, :log, :cos, :sin, :tan, :cosh, :sinh, :tanh, :acos, :asin, :atan, :acosh, :asinh, :atanh, :sqrt)
  @eval begin
    function _pullback(cx::AContext,
                       f::typeof($func),
                       A::LinearAlgebra.RealHermSymComplexHerm)
      return _pullback(cx, A -> _apply_series_func(f, A), A)
    end
  end
end

@adjoint function sincos(A::LinearAlgebra.RealHermSymComplexHerm)
  n = LinearAlgebra.checksquare(A)
  λ, U = eigen(A)
  sλ, cλ = Buffer(λ), Buffer(λ)
  for i in Base.OneTo(n)
    @inbounds sλ[i], cλ[i] = sincos(λ[i])
  end
  sinλ, cosλ = copy(sλ), copy(cλ)
  sinA, cosA = U * Diagonal(sinλ) * U', U * Diagonal(cosλ) * U'
  Ω, processback = Zygote.pullback(sinA, cosA) do s,c
    return (_process_series_matrix(sin, s, A, λ),
            _process_series_matrix(cos, c, A, λ))
  end
  return Ω, function (Ω̄)
    s̄inA, c̄osA = processback(Ω̄)
    s̄inΛ, c̄osΛ = U' * s̄inA * U, U' * c̄osA * U
    PS = _pairdiffquotmat(sin, n, λ, sinλ, cosλ, -sinλ)
    PC = _pairdiffquotmat(cos, n, λ, cosλ, -sinλ, -cosλ)
    Ā = U * (PS .* s̄inΛ .+ PC .* c̄osΛ) * U'
    return (Ā,)
  end
end

Zygote.@adjoint function LinearAlgebra.tr(x::AbstractMatrix)
  # x is a squre matrix checked by tr,
  # so we could just use Eye(size(x, 1))
  # to create a Diagonal
  tr(x), function (Δ::Number)
    (Diagonal(Fill(Δ, (size(x, 1), ))), )
  end
end

# Various sensitivities for `literal_getproperty`, depending on the 2nd argument.
@adjoint function literal_getproperty(C::Cholesky, ::Val{:uplo})
  return literal_getproperty(C, Val(:uplo)), function(Δ)
    return ((uplo=nothing, info=nothing, factors=nothing),)
  end
end
@adjoint function literal_getproperty(C::Cholesky, ::Val{:info})
  return literal_getproperty(C, Val(:info)), function(Δ)
    return ((uplo=nothing, info=nothing, factors=nothing),)
  end
end
@adjoint function literal_getproperty(C::Cholesky, ::Val{:U})
  return literal_getproperty(C, Val(:U)), function(Δ)
    Δ_factors = C.uplo == 'U' ? UpperTriangular(Δ) : LowerTriangular(copy(Δ'))
    return ((uplo=nothing, info=nothing, factors=Δ_factors),)
  end
end
@adjoint function literal_getproperty(C::Cholesky, ::Val{:L})
  return literal_getproperty(C, Val(:L)), function(Δ)
    Δ_factors = C.uplo == 'L' ? LowerTriangular(Δ) : UpperTriangular(copy(Δ'))
    return ((uplo=nothing, info=nothing, factors=Δ_factors),)
  end
end

@adjoint function logdet(C::Cholesky)
  return logdet(C), function(Δ)
    return ((uplo=nothing, info=nothing, factors=Diagonal(2 .* Δ ./ diag(C.factors))),)
  end
end

@adjoint function Matrix(::UniformScaling, i::Integer, j::Integer)
  return Matrix(I, i, j), Δ -> ((λ=tr(Δ),), nothing, nothing)
end
@adjoint function Matrix(::UniformScaling, ij::NTuple{2, Integer})
  return Matrix(I, ij), Δ -> ((λ=tr(Δ),), nothing)
end
@adjoint function Matrix{T}(::UniformScaling, i::Integer, j::Integer) where {T}
  return Matrix{T}(I, i, j), Δ -> ((λ=tr(Δ),), nothing, nothing)
end
@adjoint function Matrix{T}(::UniformScaling, ij::NTuple{2, Integer}) where {T}
  return Matrix{T}(I, ij), Δ -> ((λ=tr(Δ),), nothing)
end
@adjoint function +(A::AbstractMatrix, S::UniformScaling)
  return A + S, Δ->(Δ, (λ=tr(Δ),))
end
@adjoint function -(S::UniformScaling, A::AbstractMatrix)
  return S - A, Δ->((λ=tr(Δ),), -Δ)
end

@adjoint +(A::AbstractArray, B::AbstractArray) = A + B, Δ->(Δ, Δ)
@adjoint -(A::AbstractArray, B::AbstractArray) = A - B, Δ->(Δ, -Δ)
@adjoint -(A::AbstractArray) = -A, Δ->(-Δ,)

# FFTW
# ===================

# FFTW functions do not work with FillArrays, which are needed
# for some functionality of Zygote. To make it work with FillArrays
# as well, overload the relevant functions
FFTW.fft(x::Fill, dims...) = FFTW.fft(collect(x), dims...)
FFTW.ifft(x::Fill, dims...) = FFTW.ifft(collect(x), dims...)

# the adjoint jacobian of an FFT with respect to its input is the reverse FFT of the
# gradient of its inputs, but with different normalization factor
@adjoint function FFTW.fft(xs)
  return FFTW.fft(xs), function(Δ)
    N = length(xs)
    return (N * FFTW.ifft(Δ),)
  end
end

@adjoint function FFTW.ifft(xs)
  return FFTW.ifft(xs), function(Δ)
    N = length(xs)
    return (1/N* FFTW.fft(Δ),)
  end
end

@adjoint function FFTW.fft(xs, dims)
  return FFTW.fft(xs, dims), function(Δ)
    # dims can be int, array or tuple,
    # convert to collection for use as index
    dims = collect(dims)
    # we need to multiply by all dimensions that we FFT over
    N = prod(collect(size(xs))[dims])
    return (N * FFTW.ifft(Δ, dims), nothing)
  end
end

@adjoint function FFTW.ifft(xs,dims)
  return FFTW.ifft(xs, dims), function(Δ)
    # dims can be int, array or tuple,
    # convert to collection for use as index
    dims = collect(dims)
    # we need to divide by all dimensions that we FFT over
    N = prod(collect(size(xs))[dims])
    return (1/N * FFTW.fft(Δ, dims),nothing)
  end
end

# FillArray functionality
# =======================

@adjoint function broadcasted(op, r::AbstractFill{<:Real})
  y, _back = Zygote.pullback(op, getindex_value(r))
  back(Δ::AbstractFill) = (nothing, Fill(_back(getindex_value(Δ))[1], size(r)))
  back(Δ::AbstractArray) = (nothing, getindex.(_back.(Δ), 1))
  return Fill(y, size(r)), back
end


text:
using QXZoo
using Test
using TestSetExtensions

@testset "QXZoo.jl" begin

@includetests ARGS

end


text:
#=
 store.jl

 A track can be in

 1. mask/interval format
 2. point format

 **mask/interval** format represents intervals, for example contained in a bed file  as consecutive 1's. The
 coordinates in between the intervals contain 0's.

 **point format** represents values associated with a particular base pair coordinate. Coordinates that
 do not contain values are given the same out of range number

 The aim of the io.jl is to parse various file formats and extract the right values and put these into the store
 in either of the above formats.

 Thus far the formats that are most relevant to deal with are:

 a) methpipe's bed format for cpgs
 chr1    3000826 +       CG      1       17
 chr1    3001006 +       CG      0.8     10
 chr1    3001017 +       CG      1       12
 chr1    3001276 +       CG      1       14
 chr1    3001628 +       CG      0.867   15

 b) methpipe's bed format for hypomethylated regions
 chr1    3670331 3672761 HYPO0   184     +
 chr1    3993551 3993737 HYPO1   5       +
 chr1    4490611 4490904 HYPO2   6       +
 chr1    4491370 4497651 HYPO3   275     +
 chr1    4571560 4572248 HYPO4   58      +

 c) tab delimited format describing CpG positions
 chr1    3000826 +
 chr1    3000827 -
 chr1    3001006 +
 chr1    3001007 -
 chr1    3001017 +

 Both inserting intervals and points

 Denote:
  a) as  methpipe_bed_point
  b) as  methpipe_bed_region
  c) as  cpg_point

 Note that methpoint_point has both floating point and integer values associated with it.

 REQUIRES API:
 store_methpipe_bed_points(value_type="levels") # can also be value_type="coverage"
 store_methpipe_bed_intervals()
 store_cpg_points()

 These functions are implemented in a store-backend.jl

 We also should aim for a general form :

 # Stores Points
 store_tab_file( coord_column_idxs=[1,2], value_column_idx=[5] )

 # Stores Intervals as mask
 #  - since value is empty and three column specification for the coordinate
 store_tab_file(coord_column_idx=[1,2,3],value_column_idx=[])

 # Stores Points as mask (1 or 0) because no value specified
 store_tab_file(coord_column_idx=[1,2],value_column_idx=[])

=#

using Compat
using HDF5,JLD
using GZip
using Lumberjack
using DataFrames
include("store-backend1.jl")

function read_sample_info(sample_info_path)
    try
        if !isfile(sample_info_path)
            Lumberjack.error("The sample_info path: $sample_info_path is not a valid file")
        end
        sample_info_df = readtable(sample_info_path, allowcomments=true)
        colnames=names(sample_info_df)
        Lumberjack.info("Column names in sample info: $colnames")
        required_colnames = [:filename,:filetype,:track_id]
        for required_colname in required_colnames
            colname_appears = any(x->x==required_colname,colnames)
            if !colname_appears
               Lumberjack.error("Required column $required_colname does not appear in the sample information")
            else
               Lumberjack.info("Found required column $required_colname in sample information")
            end
        end
        return sample_info_df
    catch e
        Lumberjack.error("$e\nError reading table: Check file exists and format correct.")
    end
end

function store_samples( genomic_store_path, sample_info_path, chr_sizes_path )

    track_df=read_sample_info(sample_info_path)
    for row=1:nrow(track_df)
        track_id     = track_df[row, :track_id]
        file_type    = track_df[row, :filetype]
        filepath    = track_df[row, :filename]
        if any(x->x == :coordshift,names(track_df))
            coord_shift  = track_df[row, :coordshift]
        else
            coord_shift=0
        end
        if any(x->x == :na_val,names(track_df))
             na_val = track_df[row, :na_val]
        else
            na_val=0.0
        end

        if any(x->x == :strandfilter,names(track_df))
            strand_filter = track_df[row, :strandfilter]
        else
            strand_filter = nothing
        end

        if !isabspath( filepath )
            file_path = abspath(filepath)
        end

        if track_exists(genomic_store_path,track_id)
            Lumberjack.info("==Found track $track_id in database: skip==")
            continue
        end

        Lumberjack.info("==Saving track $track_id to database==")

        if file_type == "methpipe_bed_levels"
            store_methpipe_bed_points(filepath,
                                      genomic_store_path,
                                      track_id ,
                                      chr_sizes_path,
                                      na_val=na_val,
                                      coord_shift=coord_shift,
                                      measurement="levels")
        elseif file_type == "methpipe_bed_coverage"
            store_methpipe_bed_points(filepath ,
                                      genomic_store_path ,
                                      track_id ,
                                      chr_sizes_path,
                                      na_val=na_val,
                                      coord_shift=coord_shift,
                                      measurement="coverage")
        elseif file_type == "cpg_point"
            store_cpg_points(filepath ,
                                   genomic_store_path ,
                                   track_id ,
                                   chr_sizes_path,
                                   strand_filter_str="-",
                                   coord_shift= coord_shift)
        elseif file_type == "methpipe_bed_interval"
           store_methpipe_bed_intervals(filepath ,
                                   genomic_store_path ,
                                   track_id ,
                                   chr_sizes_path,
                                   coord_shift= coord_shift)
        elseif file_type == "bedgraph_cpg"
            store_bedgraph_cpg(filepath ,
                                   genomic_store_path ,
                                   track_id ,
                                   chr_sizes_path,
                                   coord_shift= coord_shift)
        else
            Lumberjack.info("Skip track - $track_id because filetype or storagetype not yet supported")
        end
    end
end


text:
using CodecXz
using TranscodingStreams
using Test

@testset "Xz Codec" begin
    codec = XzCompressor()
    @test codec isa XzCompressor
    @test occursin(r"^(CodecXz\.)?XzCompressor\(level=\d, check=\d+\)$", sprint(show, codec))
    @test CodecXz.initialize(codec) === nothing
    @test CodecXz.finalize(codec) === nothing

    codec = XzDecompressor()
    @test codec isa XzDecompressor
    @test occursin(r"^(CodecXz\.)?XzDecompressor\(memlimit=\d+, flags=\d+\)$", sprint(show, codec))
    @test CodecXz.initialize(codec) === nothing
    @test CodecXz.finalize(codec) === nothing

    # Generated by `lzma.compress(b"foo")` on CPython 3.5.2.
    data = Vector(b"\xfd7zXZ\x00\x00\x04\xe6\xd6\xb4F\x02\x00!\x01\x16\x00\x00\x00t/\xe5\xa3\x01\x00\x02foo\x00\x00X\x15\xa9{,\xe6,\x98\x00\x01\x1b\x03\x0b/\xb9\x10\x1f\xb6\xf3}\x01\x00\x00\x00\x00\x04YZ")
    @test read(XzDecompressorStream(IOBuffer(data))) == b"foo"
    @test read(XzDecompressorStream(IOBuffer(vcat(data, data)))) == b"foofoo"
    # corrupt data
    data[[1,3,5]] = b"bug"
    @test_throws ErrorException read(XzDecompressorStream(IOBuffer(data)))

    @test XzCompressorStream <: TranscodingStreams.TranscodingStream
    @test XzDecompressorStream <: TranscodingStreams.TranscodingStream

    TranscodingStreams.test_roundtrip_read(XzCompressorStream, XzDecompressorStream)
    TranscodingStreams.test_roundtrip_write(XzCompressorStream, XzDecompressorStream)
    TranscodingStreams.test_roundtrip_lines(XzCompressorStream, XzDecompressorStream)
    TranscodingStreams.test_roundtrip_transcode(XzCompressor, XzDecompressor)

    @test_throws ArgumentError XzCompressor(level=10)
    @test_throws ArgumentError XzDecompressor(memlimit=0)
end


text:
using Plots

using simulator.fresnelltools
include("../../src/utilities.jl")
using .utilities

n_normal = 1.4
n_spechial = 1

grating_width = 10e-9
gap_width = 400e-9

frequencies = 4e14:1e11:7e14
angles = 0:1e-2:(π/2)

layers = 100

gratings(ν, θ) = Grating(n_spechial, n_normal, grating_width, gap_width, layers, 2*π*ν/c_0, θ)
get_r_te(ν, θ) = abs((gratings(ν, θ)[1] * [1, 0])[2])^2
get_r_tm(ν, θ) = abs((gratings(ν, θ)[2] * [1, 0])[2])^2

function get_data(angles, frequencies)
    te_data = Matrix{Real}(undef, length(angles), length(frequencies))
    tm_data = Matrix{Real}(undef, length(angles), length(frequencies))
    for (i, θ) in enumerate(angles)
        for (j, ν) in enumerate(frequencies)
            te_data[i, j] = get_r_te(ν, θ)
            tm_data[i, j] = get_r_tm(ν, θ)
        end
    end
    te_data, tm_data
end

(te_data, tm_data) = get_data(angles, frequencies)
heatmap(angles, frequencies, te_data, xlabel="angles", ylabel="frequencies", title="TE modes")
heatmap(angles, frequencies, tm_data, xlabel="angles", ylabel="frequencies", title="TM modes")


text:
@rule typeof(*)(:in, Marginalisation) (m_out::PointMass, m_A::PointMass) = PointMass(mean(m_in1) * mean(m_in2))

@rule typeof(*)(:in, Marginalisation) (m_out::GammaDistributionsFamily, m_A::PointMass{ <: Real }) = begin
    return GammaShapeRate(shape(m_out), rate(m_out) * mean(m_A))
end

# if A is a matrix, then the result is multivariate
@rule typeof(*)(:in, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_A::PointMass{ <: AbstractMatrix }, meta::AbstractCorrection) = begin
    A = mean(m_A)
    ξ_out, W_out = weightedmean_precision(m_out)
    W = correction!(meta, A' * W_out * A)
    return MvNormalWeightedMeanPrecision(A' * ξ_out, W)
end

# if A is a vector, then the result is univariate
# this rule links to the special case (AbstractVector * Univariate) for forward (:out) rule
@rule typeof(*)(:in, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_A::PointMass{ <: AbstractVector }, meta::AbstractCorrection) = begin
    A = mean(m_A)
    ξ_out, W_out = weightedmean_precision(m_out)
    W = correction!(meta, dot(A, W_out, A))
    return NormalWeightedMeanPrecision(dot(A, ξ_out), W)
end

# if A is a scalar, then the input is either univariate or multivariate
@rule typeof(*)(:in, Marginalisation) (m_out::F, m_A::PointMass{ <: Real }, meta::AbstractCorrection) where { F <: NormalDistributionsFamily } = begin
    A = mean(m_A)
    ξ_out, W_out = weightedmean_precision(m_out)
    W = correction!(meta, A^2 * W_out)
    return convert(promote_variate_type(F, NormalWeightedMeanPrecision), A * ξ_out, W)
end

# specialized versions for mean-covariance parameterization
@rule typeof(*)(:in, Marginalisation) (m_out::MvNormalMeanCovariance, m_A::PointMass{ <: AbstractMatrix }, meta::AbstractCorrection) = begin
    A = mean(m_A)
    μ_out, Σ_out = mean_cov(m_out)

    z = fastcholesky(Σ_out)
    tmp = A' / z
    W = correction!(meta, tmp * A)

    return MvNormalWeightedMeanPrecision(tmp * μ_out, W)
end

@rule typeof(*)(:in, Marginalisation) (m_out::MvNormalMeanCovariance, m_A::PointMass{ <: AbstractVector }, meta::AbstractCorrection) = begin
    A = mean(m_A)
    μ_out, Σ_out = mean_cov(m_out)

    z = fastcholesky(Σ_out)
    tmp = A' / z
    W = correction!(meta, dot(tmp, A))

    return NormalWeightedMeanPrecision(dot(tmp, μ_out), W)
end