text:
module OutlierDetection
    using OutlierDetectionInterface
    const OD = OutlierDetectionInterface

    import MLJBase
    const MLJ = MLJBase

    # re-export from OutlierDetectionInterface
    export CLASS_NORMAL, CLASS_OUTLIER, DEFAULT_THRESHOLD

    # combination.jl
    export combine_mean,
           combine_median,
           combine_max

    # normalization.jl
    export scale_minmax,
           scale_unify

    # classification.jl
    export classify_quantile

    # helpers.jl
    export n_normal,
           n_outlier,
           outlier_fraction

    # mlj_helpers.jl
    export augmented_transform,
           to_categorical,
           to_univariate_finite,
           from_categorical,
           from_univariate_finite

    # mlj_wrappers.jl
    export ProbabilisticDetector,
           DeterministicDetector,
           CompositeDetector

    # mlj_transformers.jl
    export ProbabilisticTransformer,
           DeterministicTransformer,
           ScoreTransformer

    # utilities
    include("normalization.jl")
    include("classification.jl")
    include("combination.jl")
    include("helpers.jl")

    # extension
    include("mlj_helpers.jl")
    include("mlj_transformers.jl")
    include("mlj_wrappers.jl")
end


text:
using DataDeps

register(DataDep(
  "teaching-assistant-evaluation",
  "http://archive.ics.uci.edu/ml/datasets/Teaching+Assistant+Evaluation",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/tae/tae.data",
  "e83ec0a1000f4d2a7d19de9de2fc8c381357304c5c8ff942779136a9e0eefca4",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=6,
        feature_indices=1:5,
        categoric_indices=1:4,
    )
  end),
))

ext: jl

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Test, LinearMaps, LinearAlgebra, BenchmarkTools

@testset "linear combinations" begin
    CS! = LinearMap{ComplexF64}(cumsum!,
                                (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10;
                                ismutating=true)
    v = rand(ComplexF64, 10)
    u = similar(v)
    b = @benchmarkable mul!($u, $CS!, $v)
    @test run(b, samples=3).allocs == 0
    n = 10
    L = sum(fill(CS!, n))
    @test mul!(u, L, v) ≈ n * cumsum(v)
    b = @benchmarkable mul!($u, $L, $v)
    @test run(b, samples=5).allocs <= 1
    for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
        @test mul!(copy(u), L, v, α, β) ≈ Matrix(L)*v*α + u*β
    end

    A = 2 * rand(ComplexF64, (10, 10)) .- 1
    B = rand(ComplexF64, size(A)...)
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(B)
    @test @inferred(LinearMaps.MulStyle(M)) === matrixstyle
    @test @inferred(LinearMaps.MulStyle(N)) === matrixstyle
    LC = @inferred M + N
    @test @inferred(LinearMaps.MulStyle(LC)) === matrixstyle
    @test @inferred(LinearMaps.MulStyle(LC + I)) === matrixstyle
    @test @inferred(LinearMaps.MulStyle(LC + 2.0*I)) === matrixstyle
    v = rand(ComplexF64, 10)
    w = similar(v)
    b = @benchmarkable mul!($w, $M, $v)
    @test run(b, samples=3).allocs == 0
    if VERSION >= v"1.3.0-alpha.115"
        b = @benchmarkable mul!($w, $LC, $v)
        @test run(b, samples=3).allocs == 0
        for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
            if testallocs
                b = @benchmarkable mul!($w, $LC, $v, $α, $β)
                @test run(b, samples=3).allocs == 0
                b = @benchmarkable mul!($w, $(I + LC), $v, $α, $β)
                @test run(b, samples=3).allocs == 0
                b = @benchmarkable mul!($w, $(LC + I), $v, $α, $β)
                @test run(b, samples=3).allocs == 0
            end
            y = rand(ComplexF64, size(v))
            @test mul!(copy(y), LC, v, α, β) ≈ Matrix(LC)*v*α + y*β
            @test mul!(copy(y), LC+I, v, α, β) ≈ Matrix(LC + I)*v*α + y*β
        end
    end
    # @test_throws ErrorException LinearMaps.LinearCombination{ComplexF64}((M, N), (1, 2, 3))
    @test @inferred size(3M + 2.0N) == size(A)
    # addition
    @test @inferred convert(Matrix, LC) == A + B
    @test @inferred convert(Matrix, LC + LC) ≈ 2A + 2B
    @test @inferred convert(Matrix, M + LC) ≈ 2A + B
    @test @inferred convert(Matrix, M + M) ≈ 2A
    # subtraction
    @test (@inferred Matrix(LC - LC)) ≈ zeros(eltype(LC), size(LC)) atol=10eps()
    @test (@inferred Matrix(LC - M)) ≈ B
    @test (@inferred Matrix(N - LC)) ≈ -A
    @test (@inferred Matrix(M - M)) ≈ zeros(size(M)) atol=10eps()
    # scalar multiplication
    @test @inferred Matrix(-M) == -A
    @test @inferred Matrix(-LC) == -A - B
    @test @inferred Matrix(3 * M) == 3 * A
    @test @inferred Matrix(M * 3) == 3 * A
    @test Matrix(3.0 * LC) ≈ Matrix(LC * 3) ≈ 3A + 3B
    @test @inferred Matrix(3 \ M) ≈ A/3
    @test @inferred Matrix(M / 3) ≈ A/3
    @test @inferred Matrix(3 \ LC) ≈ (A + B) / 3
    @test @inferred Matrix(LC / 3) ≈ (A + B) / 3
    @test @inferred Array(3 * M' - CS!) == 3 * A' - Array(CS!)
    @test (3M - 1im*CS!)' == (3M + ((-1im)*CS!))' == M'*3 + CS!'*1im
    @test @inferred Array(M + A) == 2 * A
    @test @inferred Array(A + M) == 2 * A
    @test @inferred Array(M - A) == 0 * A
    @test Array(A - M) == 0 * A
end


text:
# Greetings! from path ~/.julia/config/startup.jl
############################################
## 1.Set parameters
############################################
# process and threads parameters
cores = 1                                                           # Int64; total number of processes. cores number should be less than "N_SR_pairs" number.
threads = 1                                                         # Int64; total number of threads. threads number should be less than N piece of data when "stack_all=false". please set to be "1" when "stack_all=true"
flag = true                                                         # true of false; print computing time for debugging purpose.
precompile = false                                                  # true or false; precompile or not.

# absolute path parameters
CCFDIR = expanduser("./CORR_fft")                              # dir where CC data are stored.
STACKDIR = expanduser("./STACK_d")                               # dir to store STACK data.

# select time-chunck folder for stack
select_time = false                                                  # true or false; whether specify time-chunck folder path.
TIMETXT = expanduser("./time_stack.txt")                            # time-chunck folder path.

# select station for stack
select_sta = "sr"                                                   # "sta" or "sr"; "sta" cross-correlate only works for selected net-sta pairs base on STATXT; "sr" cross-correlate only works for specified source and receiver based on SRTXT.
STATXT = expanduser("./sta_stack.txt")                              # station infomation including: network-name and station-name.
SRTXT = expanduser("./sr_stack.txt")                                # source and receiver information including: source-1 receiver-1 r-2 r-3...\nsource-2 r-1 r-2 r-3...
channel_regular = "*.CHN.*.CHN.*"                                   # channel(regular expression) of CCF-data; we will use 'glob(comp,CCFDIR)' to select data in 'CCFDIR' path.
comp = "NN"                                                         # component name corresponds to "channel_regular"; different from NoiseCC.CC, we will use 'load_corr(file,comp)' to read CCFData in 'CCFDIR' path, and "comp" is the group name of jld2 file(HDF5 in julia).

# stack time interval
stack_all = true                                                    # true or false; "true" for stack corr-data over all the selected time chunck.
time_interval = 1 * 24 * 60 * 60                                    # Int64 (sec); stack every "time_interval" time, and it works when stack_all = false || 中文: 每间隔多少秒叠加一次，要根据corr-data文件之间的时间间隔决定，例如corr-data文件是6个小时一个，那么time_interval=12*60*60，表示每12个小时叠加一次，即需要2个corr文件叠加在一起。该参数工作时，先搜索指定台站对的全部时间的文件，找到其中最小的开始时间与最大的结束时间，利用range函数生成间隔为""的叠加时窗，叠加每个时窗内的corr-data文件，叠加后的文件名称为该时窗的起始与终止时刻。
min_stack_chunck = 1                                                # Int64; each stack window needs at least "min_stack_chunck" to be able to stack, and it works for both stack_all = true and stack_all = false. || 中文:每个时间窗内叠加，至少需要 "min_stack_chunck" 个corr-data文件才会叠加,我们在读取corr-data文件后，会进行质量评估，剔除异常文件，导致corr-data文件数量减少。

# stack method
substack_method = "pws"                                             # "mean", "pws", and "robust"; substack method corresponding with NoiseCC.CC process, but can only choose one type.
stack_method = "mean"                                               # "mean", "pws", and "robust"; stack method, but can only choose one type.
auto_stack = true                                                   # true or false; perform auto-correlation or not.
cross_stack = true                                                  # true or false; perform cross-correlation or not.

# filter
select_filter = true                                                # true or false; whether filter before stacking.
freqmin = 0.5                                                       # Float64;
freqmax = 20.0                                                      # Float64;
corners = 4                                                         # Int64; number of corners in Butterworth filter.
zerophase = false                                                   # true or false; if true, apply filter twice for no phase shifting.

# criteria for corr-data selection
select_corr = true                                                  # true or false; whether to perform quality selecting based on median value of amp of corr-data.
median_high = 10.0                                                  # Float64; max median value (default value).
median_low = 0.0                                                    # Float64; min median value (default value).

############################################
# 2.run NoiseCC.FFT_CC
############################################
# a.using and import
# using Distributed
# NoiseCC_path = "/Users/yf/3.Project/2-1.Velocity_Change/SeisNoise_Code_YF/code_pmap_thread/"
# push!(LOAD_PATH, NoiseCC_path)
# @everywhere using NoiseCC # or @everywhere import NoiseCC

# b.include
using Distributed
addprocs(cores, topology = :master_worker, exeflags = "-t $(threads)")
@everywhere include(
    "/Users/yf/3.Project/2-1.Velocity_Change/SeisNoise_Code_YF/code_pmap_thread/NoiseCC.jl",
)

NoiseCC.STACK(
    cores,
    threads,
    flag,
    precompile,
    CCFDIR,
    STACKDIR,
    select_time,
    TIMETXT,
    select_sta,
    STATXT,
    SRTXT,
    channel_regular,
    comp,
    stack_all,
    time_interval,
    min_stack_chunck,
    substack_method,
    stack_method,
    auto_stack,
    cross_stack,
    select_filter,
    freqmin,
    freqmax,
    corners,
    zerophase,
    select_corr,
    median_high,
    median_low,
)

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


text:
## Some auxilliary functions, used for continuous wavelet analysis
## Author Bernd Blasius (bernd.blasius@gmail.com)

function deal_missing(x,t)
  # remove missing data (in csv file this is coded as NaN)
  n = length(x)
  x1 = zeros(n)
  t1 = zeros(n)

  j=0
  for i = 1:n
    if !isnan(x[i])
       j += 1
       x1[j] = x[i]
       t1[j] = t[i]
    end
  end
  x1[1:j], t1[1:j]
end

function sample_regular(x,t,dt,keep_positive)
# sample time series x,t at regular sampling rate dt and interpolate missing values
# needs julia package Dierckx

    # eliminate missing values, this yields new time vectors for every state
    xx,tx = deal_missing(x,t)

    # interpolation to a regular one-day sampling interval
    t_end = trunc(t[end])
    tt = 0:dt:t_end
    xr = evaluate(Spline1D(tx,xx),tt)   # regularly spaced time series

    # frequently we want to restrict time series to positive values
    # this may be violiated by the spline fitting which artificially may introduce
    # negative values
    # this is a simple hack to remove this, but it works (at least in the case I was interested in)
    if keep_positive
      xr = minimum(xr) < 0 ? xr .- minimum(xr) : xr
    end

    return tt,xr
end

function scramble(x,n)
# sample n random values from time series x
# to create surrogate data
    len = length(x)
    x[rand(1:len,n)]
end

function smooth(x,dt,nvoice,scale,sw=0.6)
# wavelet scale-time smoothing
# this function needs julia package DSP.jl
    xf = similar(x)  # smoothed data
    nscales, ntimes = size(x)

    # smooth in time, by filtering with a Gaussian
    for i=1:nscales
      s = scale[i]
      b = [exp(-0.5*(j*dt/s)^2)/(s*2*pi)  for j=1:ntimes]
      xf[i,:] = filtfilt(b,x[i,:])
    end

    # smooth in scale according to Si et al (2008)
    # window sw=0.6 corresponds to decorrelation length of Morlet wavelet
    sw_tot = trunc(Int,sw*nvoice)
    b = (1/sw_tot)*ones(sw_tot)     # flat (rectangular) window -> moving average
    for i=1:ntimes
      xf[:,i] = filtfilt(b,xf[:,i])
    end
    xf
 end

 function preprocess(x)
 # do some data preprocessing with a bandpass filter
    responsetype = Bandpass(1 ./ 12, 0.49; fs=1)
    designmethod = Butterworth(4)
    x = filt(digitalfilter(responsetype, designmethod), x)
    x
 end

 function preprocess1(t,x,dt)
    # bandpass filter in frequency domain (implemented by hand)
    sh = 500       #   sharpnes of filter function
    # minimum and maximum frequency cut-offs of the bandpass filter
    f1 = 1/16
    f2 = 1/5

    fmax = 1/(2*dt)           # highest possible frequency
    fmin = 1/(t[end] - t[1])  # lowest possible frequency
    freq = LinRange(fmin,2*fmax,length(t))
    # Fermi function
    fil_part1 = 1 ./ (exp.(sh*(freq .- f2)) .+ 1) .* 1. ./ (exp.(sh*(f1 .- freq)) .+ 1)
    fil_part2 = fil_part1[length(fil_part1):-1:1]  # symmetrical counter part
    fil_freq = max.(fil_part1,fil_part2)           # merge the parts

    ft = fft(x)          # Fourier transformg
    phase = angle.(ft)   # phase
    ps = abs.(ft)        # power spectrum

    #   filter data
    ps_fil = ps .* fil_freq     # filter power spectra
    # recalculate fourier
    ft_fil = cos.(phase) .* ps_fil .+       #real part
              sin.(phase) .* ps_fil*im      #imaginary part
    fil = real.(ifft(ft_fil))            # retransformation
 end

# #########################################################
# some functions for circular Statistics

function unwrap(v,vmin,vmax)
  # unwrap phase angle (phase in multiples of pi)
    if v < vmin
        v += 2.0
    elseif v > vmax
        v -= 2.0
    end
  v
end

function circstats(ph)
    # ph is a vector if phases in units of 2 pi
    n = length(ph)

    c = sum(cos.(ph)) / n
    s = sum(sin.(ph)) / n

    R = sqrt(c*c + s*s)       # order parameter
    var = 1.0 - R             # cicurlar variance
    theta = atan(s,c)         # mean direction
    std = sqrt(-2. * log(R))  # circular standard deviation
    theta, var, std
end

# #########################################################
# function to indentify strongest co-oscillating components

function get_maxInd(wco, phs, n, ind1, ind2)
    # find maximal wavelet coherence at each time instance
    # search only within a range of indices (ind1, ind2) - corresponding to large power

    wco_max_ind = zeros(Int,n)   # vector of indices of maximal coherence
    wco_max  = zeros(n)          # vector of maximal
    angle_max = zeros(n)         # vector of correspondingphase differences
    for i=1:n
        ind = argmax(wco[ind1:ind2,i]) + ind1 - 1
        wco_max_ind[i] = ind
        wco_max[i]   = wco[ind,i]
        angle_max[i] = phs[ind,i] / pi   # we always measure phase in multiples of pi
    end
    angle_max = unwrap.(angle_max, -0.5, 1.5)   # restrict phases to range -0.5 .. 1.5

    wco_max_ind, angle_max, wco_max
end

function phasehist(angles, indices, delta, width)
    # calculate histogram of phase angles (given in multiples of pi)
    delta = 0.01
    angle_hist = fit(Histogram, angles[indices],closed=:right, -0.5:delta:1.5)
    h_angle = angle_hist.weights
    h_bins = pi* (angle_hist.edges[1][2:end])

    #smoothing of histogram
    filt = (1/width)*ones(width)   # filter window
    ll = length(h_angle)
    h_angle1 = [h_angle; h_angle ; h_angle]  # wrap the histogram to consider periodic boundaries
    h_angle1 = filtfilt(filt,h_angle1)
    h_angle = h_angle1[ll+1:2*ll]

    # normalization of histogram
    #h_angle = h_angle / sum(h_angle)    # this would ormalize the histogram
    h_angle = h_angle / maximum(h_angle) # here instead we want the max to be one

    h_angle, h_bins
end

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
module FvCFD

using Printf
include("dataStructures.jl")
include("mesh.jl")
include("constitutiveRelations.jl")
include("vectorFunctions.jl")
include("timeDiscretizations.jl")
include("output.jl")
include("boundaryConditions.jl")
include("numerics.jl")
include("JST.jl")

__precompile__()

# From this file
export solve, initializeUniformSolution3D
# From boundaryConditions.jl
export supersonicInletBoundary, subsonicInletBoundary, pressureOutletBoundary, zeroGradientBoundary, wallBoundary, symmetryBoundary, emptyBoundary
# From dataStructures.jl
export Mesh, SolutionState, Fluid
# From timeDiscretizations.jl
export forwardEuler, RK2_Mid, ShuOsher, RK4, LTSEuler
# From mesh.jl
export OpenFOAMMesh

mutable struct SolverStatus
    currentTime::Float64
    nTimeSteps::Int64
    nextOutputTime::Float64
    endTime::Float64
end

######################### Initialization ###########################
# Returns cellPrimitives matrix for uniform solution
function initializeUniformSolution3D(mesh, P, T, Ux, Uy=0, Uz=0, nDims=3)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    nVars = 2 + nDims
    cellPrimitives = Matrix{Float64}(undef, nCells, nVars)
    initialValues = [ P, T, Ux, Uy, Uz ][1:nVars]

    for c in 1:nCells
        cellPrimitives[c, :] = initialValues
    end

    return cellPrimitives
end

######################### Taking and adjusting time step sizes #######################
# Calculates CFL at each cell. Expects sln.cellState, sln.cellPrimitives and sln.faceFluxes to be up to date
function CFL!(CFL, mesh::Mesh, sln::SolutionState, fluid::Fluid, dt=1)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    fill!(CFL, 0.0)
    faceRhoT = linInterp_3D(mesh, hcat(sln.cellState[:,1], sln.cellPrimitives[:,2]))
    for f in 1:nFaces
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]

        faceRho = faceRhoT[f, 1]
        faceT = faceRhoT[f, 2]

        # Use cell center values on the boundary
        if neighbourCell == -1
            faceRho = sln.cellState[ownerCell, 1]
            faceT = sln.cellPrimitives[ownerCell, 2]
        end

        @views faceVel = sln.faceFluxes[f, 1:3] ./ faceRho
        @views flux = abs(dot(faceVel, mesh.fAVecs[f]))*dt

        if faceT <= 0.0
            position = mesh.fCenters[f]
            println("Warning: Negative temperature at face $f: $position")
        end

        a = sqrt(fluid.gamma * fluid.R * faceT)
        flux += mag(mesh.fAVecs[f])*a*dt

        CFL[ownerCell] += flux
        if neighbourCell > -1
            CFL[neighbourCell] += flux
        end
    end

    CFL ./= (2 .* mesh.cVols)
end

function populateSolution(cellPrimitives, nCells, nFaces, fluid::Fluid, nDims=3)
    # Each dimension adds one momentum equation
    nConservedVars = 2+nDims
    # Each dimension adds a flux for each conserved quantity
    nFluxes = nConservedVars*nDims

    # rho, xMom, total energy from P, T, Ux, Uy, Uz
    cellState = encodePrimitives3D(cellPrimitives, fluid)
    cellFluxes = zeros(nCells, nFluxes)
    fluxResiduals = zeros(nCells, nConservedVars)
    faceFluxes = zeros(nFaces, nFluxes)

    # Initialize solution state
    sln = SolutionState(cellState, cellFluxes, cellPrimitives, fluxResiduals, faceFluxes)

    # Calculates cell fluxes, primitives from cell state
    decodeSolution_3D(sln, fluid)

    return sln
end

function restrictTimeStep(status, desiredDt)
    maxStep = min(status.endTime-status.currentTime, status.nextOutputTime-status.currentTime)

    if desiredDt > maxStep
        return true, maxStep
    else
        return false, desiredDt
    end
end

function adjustTimeStep_LTS(targetCFL, dt, status::SolverStatus)
    CFL = 1.0
    if status.nTimeSteps < 10
        # Ramp up CFL linearly in the first ten time steps to reach the target CFL
        CFL = (status.nTimeSteps+1) * targetCFL / 10
    else
        CFL = targetCFL
    end

    writeOutputThisIteration, CFL = restrictTimeStep(status, CFL)

    # Store the target CFL for the present time step in the first element of the dt vector.
    # The time discretization function will determine the actual local time step based on this target CFL
    dt[1] = CFL # TODO: Cleaner way to pass this information

    return writeOutputThisIteration, dt, CFL
end

function adjustTimeStep(maxCFL, targetCFL, dt, status)
    # If CFL too high, attempt to preserve stability by cutting timestep size in half
    if maxCFL > targetCFL*1.01
        dt *= targetCFL/(2*maxCFL)
    # Otherwise slowly approach target CFL
    else
        dt *= ((targetCFL/maxCFL - 1)/10+1)
    end

    writeOutputThisIteration, dt = restrictTimeStep(status, dt)
    return writeOutputThisIteration, dt, maxCFL
end

function advanceStatus!(status::SolverStatus, dt, CFL, timeIntegrationFn, silent)
    status.nTimeSteps += 1

    if timeIntegrationFn == LTSEuler
        status.currentTime += CFL
    else
        status.currentTime += dt
    end

    if !silent
        @printf("Timestep: %5.0f, simTime: %9.4g, Max CFL: %9.4g \n", status.nTimeSteps, status.currentTime, CFL)
    end
end

defaultAirFluid = Fluid(1005, 287.05, 1.4)

######################### Main #######################
#=
    This is where the magic happens!
    Do CFD!

    Arguments:
        mesh: see dataStructuresDefinitions.md / dataStructures.jl
        meshPath:           (string) path to folder where mesh is stored
        cellPrimitives:     initial cell-center conditions, see dataStructuresDefinitions.md / dataStructures.jl
        boundaryConditions: Array of alternating boundary condition function references and associated boundary condition parameters:
            Ex: [ emptyBoundary, [], supersonicInletBoundary, [P, T, Ux, Uy, Uz], zeroGradientBoundary, [], symmetryBoundary, [], zeroGradientBoundary, [], wallBoundary, [] ]
            Order of BC's must match the order of boundaries defined in the mesh's 'boundaries' file
        timeIntegrationFn:  The desired time integration function from timeDiscretizations.jl (should update sln.cellState, sln.cellPrimitives, and sln.cellFluxes based on their old values and sln.cellResiduals (obtained from fluxFunction))
        fluxFunction:       The desired function to calculate sln.cellResiduals from sln.cellState, sln.cellPrimitives, and sln.cellFluxes
        initDt:             Initial time step (s)
        endTime:            Simulation End Time (s). Start Time = 0
        outputInterval:     Writes solution/restart files every outputInterval seconds of simulation time
        fluid:              Fluid definition of type Fluid (see dataStructuresDefinitions.md / dataStructures.jl)
        targetCFL:          Target maximum CFL in the computation domain (time step will be adjusted based on this value)
        gamma/R/Cp:         Fluid properties
        silent:             Controls whether progress is written to console (can slow down simulation slightly for very simple cases)
        restart:            Controls whether restarting from a restart file
        createRestartFile:  Controls whether to write restart files. These are overwritten every time they are outputted.
        createVTKOutput:    Controls whether to write .vtk output. These are not overwritten every time they are outputted.
        restartFiles:       (string) path to restart file to read/write from

    Returns:
        sln.cellPrimitives at end of simulation

    Outputs:
        restart file (if specified)
        .vtk files (is specified)
=#
function solve(mesh::Mesh, meshPath, cellPrimitives::Matrix{Float64}, boundaryConditions, timeIntegrationFn=forwardEuler,
        fluxFunction=unstructured_JSTFlux; initDt=0.001, endTime=0.14267, outputInterval=0.01, targetCFL=0.2, fluid=defaultAirFluid,
        silent=true, restart=false, createRestartFile=true, createVTKOutput=true, restartFile="FvCFDRestart.txt")

    if !silent
        println("Initializing Simulation")
    end

    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    if !silent
        println("Mesh: $meshPath")
        println("Cells: $nCells")
        println("Faces: $nFaces")
        println("Boundaries: $nBoundaries")
    end

    if restart
        if !silent
            println("Reading restart file: $restartFile")
        end

        cellPrimitives = readRestartFile(restartFile)

        # Check that restart file has the same number of cells as the mesh
        nCellsRestart = size(cellPrimitives, 1)
        if nCellsRestart != nCells
            throw(ArgumentError("Number of cells in restart file ($nCellsRestart) does not match the present mesh ($nCells). Please provide a matching mesh and restart file."))
        end
    end

    sln = populateSolution(cellPrimitives, nCells, nFaces, fluid, 3)

    dt = initDt
    if timeIntegrationFn==LTSEuler
        # If using local time stepping, the time step will be different for each cell
        dt = zeros(nCells)
    end

    status = SolverStatus(0, 0, outputInterval, endTime)
    writeOutputThisIteration = false
    CFLvec = zeros(nCells)

    if !silent
        println("Starting iterations")
    end

    ### Main Loop ###
    while status.currentTime < status.endTime
        ############## Adjust timestep size #############
        if timeIntegrationFn == LTSEuler # LTS = Local time stepping
            writeOutputThisIteration, dt, CFL = adjustTimeStep_LTS(targetCFL, dt, status)
        else
            CFL!(CFLvec, mesh, sln, fluid, dt)
            writeOutputThisIteration, dt, CFL = adjustTimeStep(maximum(CFLvec), targetCFL, dt, status)
        end

        ############## Take a timestep #############
        sln = timeIntegrationFn(mesh, fluxFunction, sln, boundaryConditions, fluid, dt)
        advanceStatus!(status, dt, CFL, timeIntegrationFn, silent)

        ############## Write Output #############
        if writeOutputThisIteration
            writeOutput(sln.cellPrimitives, restartFile, meshPath, createRestartFile, createVTKOutput)
            status.nextOutputTime += outputInterval
        end
    end

    # Always create output upon exit
    writeOutput(sln.cellPrimitives, restartFile, meshPath, createRestartFile, createVTKOutput)

    # Return current cell-center properties
    return sln.cellPrimitives
end

end

max_stars_repo_name: henrystoldt/fvCFD


text:
include("./FieldFactory/rayclassgrp.jl")
include("./FieldFactory/ab_exts.jl")


text:
using Test
using DotNET
import DotNET: unbox, arrayof

@testset "Type Loader" begin
    @test T"System.Int64".Name == "Int64"
    @test T"System.String".Name == "String"
    @test T"System.Console, mscorlib".Name == "Console"
    @test_throws CLRException T"BadType.KUUPXXCL"
end # testset

@testset "Marshaller" begin
    uturn(x) = @test (unbox(convert(CLRObject, x)) === x)
    true
    uturn(Int8(42))
    uturn(UInt8(42))
    uturn(Int16(42))
    uturn(UInt16(42))
    uturn(Int32(42))
    uturn(UInt32(42))
    uturn(Int64(42))
    uturn(UInt64(42))
    uturn(Float32(42))
    uturn(Float64(42))
    uturn('A')
    uturn("Hello, World!")
    array = convert(CLRObject, [1, 2, 3])
    @test collect(array) == [1, 2, 3]
    for (i, x) in enumerate(array)
        @test x == i
    end
end

@testset "Member Invocation" begin
    ArrayList = T"System.Collections.ArrayList, mscorlib"
    li = ArrayList.new()
    li.Add(42)
    @test li.Count == 1
    ListT = T"System.Collections.Generic.List`1, mscorlib"
    li = ListT.new[T"System.Int64"]()
    li.Add(42)
    @test li.Count == 1
    @test isclrtype(T"System.Array".Empty[T"System.Int64"](), T"System.Int64[]")
    WeakReference = T"System.WeakReference, mscorlib"
    ref = WeakReference.new(CLRObject(0))
    ref.Target = WeakReference
    @test T"System.Object, mscorlib".Equals(ref.Target, WeakReference)

    # varargs
    @test T"System.String".Format("{0}{1}{2}", "i", 18, "n") == "i18n"

    # field access
    let t = T"System.ValueTuple`1".new[T"System.Int32"](Int32(1))
        @test t.Item1 == Int32(1)
        t.Item1 = Int32(2)
        @test t.Item1 == Int32(2)
    end

    # out parameters
    let
        tryparse(s, i) = T"System.Int64".TryParse(s, i)
        @test tryparse("42", null)
        @test !tryparse("abc", null)
        r = Ref{Int64}()
        @test !tryparse("abc", r)
        tryparse("123", r)
        @test r[] == 123
        r = Ref{CLRObject}()
        tryparse("def", r)
        @test unbox(r[]) == 0
        tryparse("456", r)
        @test unbox(r[]) == 456
        r = Ref{CLRObject}(1)
        tryparse("789", r)
        @test unbox(r[]) == 789
    end
end

@testset "Generics" begin
    List = T"System.Collections.Generic.List`1"
    li = List.new[T"System.Int64"]()
    li.Add(Int64(1))
    @test eltype(collect(li)) == Int64
end

@testset "Arrays and Iteration" begin
    @test eltype(convert(CLRObject, Int64[1])) == Int64
    jlarr = reshape(1:24, (2, 3, 4))
    arr = convert(CLRObject, jlarr)
    @test axes(arr) == (2, 3, 4)
    @test collect(arr) == jlarr
    ch = Channel() do ch
        e = arr.GetEnumerator()
        while e.MoveNext()
            put!(ch, e.Current)
        end
    end
    @test collect(ch) == jlarr[:]
end


text:
using Flux
using Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated

include("C:\\Dropbox\\julialangfiles\\ANNBN\\ANNBN.jl")
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, MLDatasets, Plots, Printf, Debugger
train_x, yy_train_all = MLDatasets.MNIST.traindata()
test_x,  yy_test_all  = MLDatasets.MNIST.testdata()
xx_train=MLDatasets.MNIST.convert2features(MLDatasets.MNIST.traintensor())'
xx_test=MLDatasets.MNIST.convert2features(test_x)'
xx_train=convert(Array{Float64,2},xx_train)
xx_test=convert(Array{Float64,2},xx_test)
yy_train_all=convert(Array{Float64,1},yy_train_all)
yy_test_all=convert(Array{Float64,1},yy_test_all)
i_train=size(xx_train,1)
i_test=size(xx_test,1)
vars=size(xx_train,2)

tol1=0.0
num_to_test=0
yy_train=copy(yy_train_all)
i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1

yy_test=copy(yy_test_all)
i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1

# neurons=Int64(floor(i_train/(vars+1)))
# neurons=100
# inds_all,n_per_part=ann_by_parts.___clustering(neurons,xx_train,200)
# x1=copy(xx_train[inds_all[n_per_part[1]+1:n_per_part[1+1]],:])
# obs1=size(x1,1)
# phi1=ANNBN.calc_phi_deriv(x1,x1,obs1,obs1,cc1,vars,rand(1:vars))
# phi_deriv=ANNBN.calc_phi(x1,x1,obs1,obs1,cc1,vars)
# aa=[phi1;phi_deriv]\[yy_train[inds_all[n_per_part[1]+1:n_per_part[1+1]]];zeros(obs1)]
# mean(abs.([phi1;phi_deriv]*aa-[yy_train[inds_all[n_per_part[1]+1:n_per_part[1+1]]];zeros(obs1)]))

# i_train=10000
include("C:\\Dropbox\\julialangfiles\\ANNBN\\ANNBN.jl")
inds_all=1:i_train;neurons=10000;items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
cc1=1.0
a_all,itmps=ANNBN.train_layer_1_rbf_deriv(neurons,vars,i_train,n_per_part,inds_all,xx_train[1:i_train,:],yy_train[1:i_train],cc1)
# phi_tr=ANNBN.calc_phi(xx_train[1:i_train,:],xx_train[1:i_train,:],i_train,i_train,cc1,vars)
layer1=Array{Float64}(undef,i_train,0)
@time for i=1:length(a_all)
    global layer1,inds_all,n_per_part,a_all,phi_tr
    # layer1=[layer1 phi_tr[:,inds_all[n_per_part[i]+1:n_per_part[i+1]]]*a_all[i]]
    layer1=[layer1 phi_tr[:,itmps[i]]*a_all[i]]
    println(i)
end
# @time a_layer1=inv([layer1 ones(i_train)]'*[layer1 ones(i_train)])*[layer1 ones(i_train)]'*yy_train[1:i_train]
@time a_layer1=[layer1 ones(i_train)]\yy_train[1:i_train]
predl1=[layer1 ones(i_train)]*a_layer1
maetr=mean(abs.(yy_train[1:i_train]-predl1))
maximum(abs.(yy_train[1:i_train]-predl1))

# phi_te=ANNBN.calc_phi(xx_train,xx_test,i_train,i_test,cc1,vars)
layer1_test=Array{Float64}(undef,i_test,0)
for i=1:length(a_all)
    global layer1_test,inds_all,n_per_part,a_all,phi_te
    # layer1_test=[layer1_test phi_te[:,inds_all[n_per_part[i]+1:n_per_part[i+1]]]*a_all[i]]
    layer1_test=[layer1_test phi_te[:,itmps[i]]*a_all[i]]
end
predl1_test=[layer1_test ones(i_test)]*a_layer1
maete=mean(abs.(yy_test-predl1_test))
maximum(abs.(yy_test-predl1_test))

i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train[1:i_train]).<0.01;maetr=100sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=100sum(i2)/length(i2)

using JLD2, JSON
@save "phis.jld2" phi_tr phi_te
@load "phis.jld2" phi_tr phi_te

##### APPROXIMATION OF DERIVATIVES
# first variable
yy_train_deriv=-3xx_train[:,1].^2+10xx_train[:,1]+5xx_train[:,2]
yy_train_deriv_deriv=-6xx_train[:,1]+10ones(i_train)
yy_test_deriv=-3xx_test[:,1].^2+10xx_test[:,1]+5xx_test[:,2]
yy_test_deriv_deriv=-6xx_test[:,1]+10ones(i_test)

predl1_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_deriv-yy_train_deriv))
predl1_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_deriv_deriv-yy_train_deriv_deriv))
predl1_test_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_test_deriv-yy_test_deriv))
predl1_test_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_test_deriv_deriv-yy_test_deriv_deriv))

# second variable
yy_train_deriv=5xx_train[:,1]
yy_train_deriv_deriv=zeros(i_train)
yy_test_deriv=5xx_test[:,1]
yy_test_deriv_deriv=zeros(i_test)

predl1_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_deriv-yy_train_deriv))
predl1_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_deriv_deriv-yy_train_deriv_deriv))

predl1_test_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_test_deriv-yy_test_deriv))
predl1_test_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_test_deriv_deriv-yy_test_deriv_deriv))


text:
using Celeste: DeterministicVI, SensitiveFloats
import DeterministicVI: calculate_G_s!, load_source_brightnesses, elbo_likelihood

using Base.Test
using Distributions

import ForwardDiff.Dual

import SampleData: gen_two_body_dataset, true_star_init

@testset "elbo" begin

@testset "calculate_G_s overwrites E_G_s, E_G2_s, and var_G_s" begin
    ea, vp, catalog = gen_two_body_dataset()
    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, ea.Sa, true, true)

    # let's write some non-trivial derivatives to the fs0m and fs1m
    elbo_likelihood(ea, vp, elbo_vars)

    sbs = load_source_brightnesses(ea, vp)

    # for this call E_G_s and E_G2_s are initialized to zero
    calculate_G_s!(vp, elbo_vars, sbs[2], 4, 2, true)
    ev_cleared = deepcopy(elbo_vars)

    # for this call E_G_s and E_G2_s have not been zeroed out
    calculate_G_s!(vp, elbo_vars, sbs[1], 2, 1, true)
    calculate_G_s!(vp, elbo_vars, sbs[2], 3, 2, true)
    calculate_G_s!(vp, elbo_vars, sbs[2], 4, 2, true)

    @test ev_cleared.E_G_s.v[] ≈ elbo_vars.E_G_s.v[]
    @test ev_cleared.E_G_s.d[] ≈ elbo_vars.E_G_s.d[]
    @test ev_cleared.E_G_s.h[] ≈ elbo_vars.E_G_s.h[]

    @test ev_cleared.E_G2_s.v[] ≈ elbo_vars.E_G2_s.v[]
    @test ev_cleared.E_G2_s.d[] ≈ elbo_vars.E_G2_s.d[]
    @test ev_cleared.E_G2_s.h[] ≈ elbo_vars.E_G2_s.h[]

    @test ev_cleared.var_G_s.v[] ≈ elbo_vars.var_G_s.v[]
    @test ev_cleared.var_G_s.d[] ≈ elbo_vars.var_G_s.d[]
    @test ev_cleared.var_G_s.h[] ≈ elbo_vars.var_G_s.h[]
end

@testset "test bvn cov" begin
    gal_axis_ratio = .7
    gal_angle = pi/5
    gal_radius_px = 2.

    manual_11 = gal_radius_px^2 * (1 + (gal_axis_ratio^2 - 1) * (sin(gal_angle))^2)
    util_11 = DeterministicVI.get_bvn_cov(gal_axis_ratio, gal_angle, gal_radius_px)[1,1]
    @test util_11 ≈ manual_11

    manual_12 = gal_radius_px^2 * (1 - gal_axis_ratio^2) * (cos(gal_angle)sin(gal_angle))
    util_12 = DeterministicVI.get_bvn_cov(gal_axis_ratio, gal_angle, gal_radius_px)[1,2]
    @test util_12 ≈ manual_12

    manual_22 = gal_radius_px^2 * (1 + (gal_axis_ratio^2 - 1) * (cos(gal_angle))^2)
    util_22 = DeterministicVI.get_bvn_cov(gal_axis_ratio, gal_angle, gal_radius_px)[2,2]
    @test util_22 ≈ manual_22
end

@testset "test active sources" begin
    # Test that the derivatives of the expected brightnesses partition in
    # active_sources.
    ea, vp, catalog = gen_two_body_dataset()

    # for subsequent tests, ensure that the first light source
    # has a radius at least as large as both the height and the
    # width of images, and that its center is within the image
    s = 1
    for n = 1:size(ea.patches, 2)
        p = ea.patches[s, n]
        @test 0.5 <= p.world_center[1] <= 20.5
        @test 0.5 <= p.world_center[2] <= 23.5

        # this patch is huge, the bottom left corner should be the
        # bottom left of the image
        @test p.bitmap_offset == [0, 0]

        # all the pixels should be active pixels
        @test size(p.active_pixel_bitmap) == size(ea.images[n].pixels)
        @test all(p.active_pixel_bitmap)
    end

    n = 5
    p = ea.patches[s, n]

    # lets make the second source have no active pixels
    s2 = 2
    for i in 1:ea.N
        fill!(ea.patches[s2, i].active_pixel_bitmap, false)
    end

    P = length(CanonicalParams)
    ea_no2 = ElboArgs(ea.images, ea.patches, [1, 2])
    elbo_no2 = DeterministicVI.elbo_likelihood(ea_no2, vp)
    @test elbo_no2.d[:, 2] ≈ zeros(P, 1)

    # lets make the second source have only a few active pixels
    # in one band, so it and source 1 don't overlap completely
    p2 = ea.patches[s2, n]
    p2.active_pixel_bitmap[10:11,10:11] = true

    # now on to a main test--active source order shouldn't matter
    ea12 = ElboArgs(ea.images, ea.patches, [1, 2])
    elbo_lik_12 = DeterministicVI.elbo_likelihood(ea12, vp)

    ea21 = ElboArgs(ea.images, ea.patches, [2, 1])
    elbo_lik_21 = DeterministicVI.elbo_likelihood(ea21, vp)
    @test elbo_lik_12.v[] ≈ elbo_lik_21.v[]
    @test elbo_lik_12.d[:,1] ≈ elbo_lik_21.d[:,2]
    @test elbo_lik_12.d[:,2] ≈ elbo_lik_21.d[:,1]

    # next main test--deriviatives for active sources don't
    ea1 = ElboArgs(ea.images, ea.patches, [1,])
    elbo_lik_1 = DeterministicVI.elbo_likelihood(ea1, vp)
    # source 1 includes all of source 2's active pixels
    @test elbo_lik_1.v[] ≈ elbo_lik_12.v[]

    ea2 = ElboArgs(ea.images, ea.patches, [2,])
    elbo_lik_2 = DeterministicVI.elbo_likelihood(ea2, vp)

    @test elbo_lik_12.d[:, 1] ≈ elbo_lik_1.d[:, 1]
    @test elbo_lik_12.d[:, 2] ≈ elbo_lik_2.d[:, 1]

    @test elbo_lik_12.h[1:P, 1:P] ≈ elbo_lik_1.h
    @test elbo_lik_12.h[(1:P) + P, (1:P) + P] ≈ elbo_lik_2.h
end

@testset "star truth is most likely" begin
    ea, vp, catalog = true_star_init()

    best = DeterministicVI.elbo_likelihood(ea, vp)

    for bad_a in [.3, .5, .9]
        vp_a = deepcopy(vp)
        vp_a[1][ids.is_star] = [ 1.0 - bad_a, bad_a ]
        bad_a_lik = DeterministicVI.elbo_likelihood(ea, vp_a)
        @test best.v[] > bad_a_lik.v[]
    end

    for h2 in -2:2
        for w2 in -2:2
            if !(h2 == 0 && w2 == 0)
                vp_mu = deepcopy(vp)
                vp_mu[1][ids.pos] += [h2 * .5, w2 * .5]
                bad_mu = DeterministicVI.elbo_likelihood(ea, vp_mu)
                @test best.v[] > bad_mu.v[]
            end
        end
    end

    for delta in [.7, .9, 1.1, 1.3]
        vp_r1 = deepcopy(vp)
        vp_r1[1][ids.flux_loc] += log(delta)
        bad_r1 = DeterministicVI.elbo_likelihood(ea, vp_r1)
        @test best.v[] > bad_r1.v[]
    end

    for b in 1:4
        for delta in [-.3, .3]
            vp_c1 = deepcopy(vp)
            vp_c1[1][ids.color_mean[b, 1]] += delta
            bad_c1 = DeterministicVI.elbo_likelihood(ea, vp_c1)
            @test best.v[] > bad_c1.v[]
        end
    end
end

@testset "galaxy truth is most likely" begin
    ea, vp, catalog = gen_sample_galaxy_dataset(perturb=false)
    vp[1][ids.is_star] = [ 0.01, .99 ]
    best = DeterministicVI.elbo_likelihood(ea, vp)

    for bad_a in [.3, .5, .9]
        vp_a = deepcopy(vp)
        vp_a[1][ids.is_star] = [ 1.0 - bad_a, bad_a ]
        bad_a = DeterministicVI.elbo_likelihood(ea, vp_a)
        @test best.v[] > bad_a.v[];
    end

    for h2 in -2:2
        for w2 in -2:2
            if !(h2 == 0 && w2 == 0)
                vp_mu = deepcopy(vp)
                vp_mu[1][ids.pos] += [h2 * .5, w2 * .5]
                bad_mu = DeterministicVI.elbo_likelihood(ea, vp_mu)
                @test best.v[] > bad_mu.v[]
            end
        end
    end

    for bad_scale in [.8, 1.2]
        vp_r1 = deepcopy(vp)
        vp_r1[1][ids.flux_loc] += 2 * log(bad_scale)
        bad_r1 = DeterministicVI.elbo_likelihood(ea, vp_r1)
        @test best.v[] > bad_r1.v[]
    end

    for n in [:gal_axis_ratio, :gal_angle, :gal_radius_px]
        for bad_scale in [.8, 1.2]
            vp_bad = deepcopy(vp)
            vp_bad[1][getfield(ids, n)] *= bad_scale
            bad_elbo = DeterministicVI.elbo_likelihood(ea, vp_bad)
            @test best.v[] > bad_elbo.v[]
        end
    end

    for b in 1:4
        for delta in [-.3, .3]
            vp_c1 = deepcopy(vp)
            vp_c1[1][ids.color_mean[b, 2]] += delta
            bad_c1 = DeterministicVI.elbo_likelihood(ea, vp_c1)
            @test best.v[] > bad_c1.v[]
        end
    end
end

@testset "manual elbo gradient matches auto diff gradient" begin
    ea, vp, catalog = gen_two_body_dataset()

    # compute the gradient (and the value) of the elbo manually
    elbo_float = DeterministicVI.elbo(ea, vp)

    # create elbo arguments of the Dual number type, that compute the gradient but not
    # the hessian (it doesn't matter that the "perturbation" are all zero, it's just
    # for testing the speed and verifying that it works)
    vp_dual = convert(VariationalParams{Dual{1, Float64}}, vp)

    for s in 1:2, i in 1:length(ids)
        vp_dual[s][i] += ForwardDiff.Dual(0, 1)
        elbo_dual = DeterministicVI.elbo(ea, vp_dual)
        vp_dual[s][i] -= ForwardDiff.Dual(0, 1)

        fwd_diff = elbo_dual.v[].partials[]
        manual_diff = elbo_float.d[i, s]
        @test manual_diff ≈ fwd_diff
    end
end

@testset "manual elbo hessian column sums match auto diff" begin
    ea, vp, catalog = gen_two_body_dataset()

    # compute the hessian of the elbo manually
    elbo_float = DeterministicVI.elbo(ea, vp)

    # create elbo arguments of the Dual number type, that compute the gradient but not
    # the hessian (it doesn't matter that the "perturbation" are all zero, it's just
    # for testing the speed and verifying that it works)
    vp_dual = convert(VariationalParams{Dual{1, Float64}}, vp)

    for s in 1:2, i in 1:length(ids)
        # we're effectively multiplying the hessian by the 1's vector here.
        vp_dual[s][i] += ForwardDiff.Dual(0, 1)
    end

    elbo_dual = DeterministicVI.elbo(ea, vp_dual)

    for s in 1:2, p in 1:length(ids)
        auto_hessian_column_sum = elbo_dual.d[p, s].partials[]
        col_id = (s - 1) * length(ids) + p
        manual_hessian_column_sum = sum(elbo_float.h[:, col_id])
        @test manual_hessian_column_sum ≈ auto_hessian_column_sum
    end
end

@testset "automatic hessian vector product matches finite diff" begin
    ea, vp, catalog = gen_two_body_dataset()

    # compute the gradient (and the value) of the elbo manually
    d0 = DeterministicVI.elbo(ea, vp).d[:]

    vp1 = deepcopy(vp)
    perturbation = 1e-5
    vp1[1] += perturbation
    vp1[2] += perturbation
    d1 = DeterministicVI.elbo(ea, vp1).d[:]

    hv_manual = (d1 - d0) / perturbation

    vp_dual = convert(VariationalParams{Dual{1, Float64}}, vp)

    for s in 1:2, i in 1:length(ids)
        vp_dual[s][i] += ForwardDiff.Dual(0, 1)
    end

    elbo_dual = DeterministicVI.elbo(ea, vp_dual)

    P = length(ids)
    hv_auto = [elbo_dual.d[i].partials[] for i in 1:(2P)]

    for i in 1:20
        @test hv_manual[i] ≈ hv_auto[i] atol=abs(0.01 * hv_auto[i])
    end
end

end


text:
# SPDX-License-Identifier: MIT

"""
    fstprint(file, isyms=\"\", osyms=\"\")

Prints out binary FSTs in simple text format.

You can provide `isyms` and/or `osyms` mapping file (e.g. Kaldi's words.txt) to map input/output symbol IDs into symbols (e.g. words).
"""
fstprint(file; isyms="", osyms="") = `fstprint --isymbols=$isyms --osymbols=$osyms $file`

"""
	draw_fst(filename, isyms=\"\", osyms=\"\")

Wrapper around FST draw command, and graphviz dot, that displays the resulting visualisation.
"""
function draw_fst(filename; isyms="", osyms="")
    img = read(pipeline(
        `fstdraw --isymbols=$isyms --osymbols=$osyms --portrait $filename`,
        `dot -Tsvg`
    ), String)
    display("image/svg+xml", img)
end

function read_wfst(SF::Type{<:Semifield}, wfstfile; isyms="", osyms="")
	fst = read(fstprint(wfstfile; isyms=isyms, osyms=osyms), String);
	fst_finals = Dict()
	fst_arcs = Dict()

	for line in split(fst, '\n')
	    arc = split(line, '\t')
	    if length(arc) < 2
	        continue
	    elseif length(arc) < 4
	        # terminal state
	        final_node = parse(Int,arc[1])
	        final_weight = parse(Float64, arc[2])
	        fst_finals[final_node] = SF(-final_weight)
	    else
	        # non-terminal arc
	        from_node = parse(Int,arc[1])
	        to_node = parse(Int,arc[2])
	        isym = parse(Int, arc[3]) # transition ID
	        osym = arc[4] # word
	        weight = if length(arc) == 5 parse(Float64, arc[5]) else 0.0 end
	        fst_arcs[from_node] = push!(get(fst_arcs, from_node, []), (to_node, isym, osym, SF(-weight)))
	    end
	end

	return fst_arcs, fst_finals
end

function load_from_wfst!(fsm::AbstractFSM{SF}, wfstfile::String; symbol_table::String = "") where SF<:Semifield
	fst_arcs, fst_finals = read_wfst(SF, wfstfile, osyms=symbol_table)
	fst2fsm = Dict()

	s = addstate!(fsm, nothing, initweight = one(SF))
	fst2fsm[0] = s
	to_visit = [0]
	visited = []

	while !isempty(to_visit)
	    from_node = pop!(to_visit)
	    push!(visited, from_node)
	    for (to_node, isym, osym, weight) in fst_arcs[from_node]
	        if !(to_node in visited) && !(to_node in to_visit)
	            push!(to_visit, to_node)
                label = if string(isym) != "0" "$isym:$osym" else nothing end
	            fst2fsm[to_node] = addstate!(fsm, label, finalweight=get(fst_finals, to_node, zero(SF)))
	        end
	        addarc!(fsm, fst2fsm[from_node], fst2fsm[to_node], weight)
	    end
	end

	return fsm
end

"""
	load_from_kaldi(wfst_file, model_file, symbol_table= \"\")

Load Kaldi recognition network from `wfst_file`.

`model_file` (e.g. path/to/final.mdl) should contain Kaldi's transition model with the mapping transition ID -> pdf ID. This is used to create `MatrixFSM`.

`symbol_table` (e.g. path/to/graph/words.txt) is output symbol mapping table.
"""
function load_from_kaldi(wfst_file::String, kaldi_model::String; symbol_table::String = "")
    SF = LogSemifield{Float64}
    fsm = VectorFSM{SF}()
    MarkovModels.load_from_wfst!(fsm, wfst_file, symbol_table=symbol_table)

    trans_mdl = open("/mnt/matylda3/ikocour/tools/kaldi/egs/tidigits/s5/exp/mono0a/final.mdl") do fd
        Kaldi.is_binary(fd)
        Kaldi.load_transition_model(fd)
    end

	isselfloop(entry, hmm_state, transition_index) = begin
	    transition_index <= length(entry[hmm_state].transitions) &&
			entry[hmm_state].transitions[transition_index].index == hmm_state
	end

	tid2pdf = Dict() # transition ID -> pdf ID mapping
	cur_transition_id = 0
	# Load tid2pdf
	for tstate in trans_mdl.tuples
	    hmm_state = tstate.hmm_state+1
	    entry = trans_mdl.topo[tstate.phone].entry
	    transitions = entry[hmm_state].transitions
	    prev_transition_id = cur_transition_id
	    cur_transition_id += length(transitions)

	    for tid in prev_transition_id+1:(cur_transition_id)
	        if isselfloop(entry, hmm_state, tid - prev_transition_id) # self-loop
	            tid2pdf[tid] = tstate.self_pdf + 1 # Kaldi index from 0
	        else # forward-step
	            tid2pdf[tid] = tstate.forward_pdf + 1 # Kaldi index from 0
	        end
	    end
	end

    fsm_no_eps = MarkovModels.remove_label(fsm, nothing)
	keyfn(label) = split(label, ":") |> first |> x -> parse(Int, x)
	return MatrixFSM(fsm_no_eps, tid2pdf, keyfn)
end


text:
module Logistic

include("LogisticTwoSpecies.jl")
export TwoSpecies

end


text:
@doc """
Dipole type.

#### Parameters

* coord_system: The coordinate system that the locations are stored in
* x,y,z: Location of dipole
* x,y,z/ori: Orientation of dipole
* color: Color of dipole for plotting
* state: State of dipol
* size: size of dipole

""" ->
type Dipole
    coord_system::AbstractString
    x::SIUnits.SIQuantity{AbstractFloat,1,0,0,0,0,0,0,0,0} # m
    y::SIUnits.SIQuantity{AbstractFloat,1,0,0,0,0,0,0,0,0} # m
    z::SIUnits.SIQuantity{AbstractFloat,1,0,0,0,0,0,0,0,0} # m
    xori::Number
    yori::Number
    zori::Number
    color::Number
    state::Number
    size::Number
end

import Base.show
function Base.show(io::IO, d::Dipole)
    @printf("Dipole with coordinates x = % 6.2f m, y = % 6.2f m, z = % 6.2f m, size = % 9.5f\n",
        float(d.x), float(d.y), float(d.z), float(d.size))
end

function Base.show(io::IO, dips::Array{Dipole})
    @printf("%d dipoles\n", length(dips))
    for d in dips
        @printf("  Dipole with coordinates x = % 6.2f m, y = % 6.2f m, z = % 6.2f m and size = % 9.5f\n",
            float(d.x), float(d.y), float(d.z), float(d.size))
    end
end


text:
using PolarimetryCalibration, Plots, LaTeXStrings, JLD, NLopt, Debug

cal = load("../data/3C286_20160229.jld")
tovector(a) = reshape(a, length(a))

pa = tovector(cal["pa"])
r = abs(pa) .< 55.0
pa = pa[r]
I = tovector(cal["i"][r])
Q = tovector(cal["q"][r])
U = tovector(cal["u"][r])
V = tovector(cal["v"][r])

PArange = minimum(pa):1.0:maximum(pa)

function cost_circ(x, grad)
    sumsq = 0.0

    for ind in 1:length(pa)
        model = müller_rec(ϵ = x[1], ϕ = x[2], ΔG = x[3], ψ = x[4],
                           α = deg2rad(45.0))* # Circular feed
#                           α = x[7])* # Circular feed
        müller_sky(pa[ind])*
#        lin_pol_source(flux = x[5], p = x[6], pa = deg2rad(34.0))
        lin_pol_source(flux = x[5], p = x[6], pa = x[7])
#        lin_pol_source(flux = 0.022, p = x[6], pa = x[7])
#        sumsq += (model[1] - I[ind])^2 +
        sumsq += (model[1] - I[ind])^2 +
        1/x[6]*(model[2] - Q[ind])^2 +
        1/x[6]*(model[3] - U[ind])^2 +
        (model[4] - V[ind])^2
#        sumsq += (model[3] - U[ind])^2 + (model[4] + Q[ind])^2
    end

    pmodel = Float64[(müller_rec(ϵ = x[1],
                                 ϕ = x[2],
                                 ΔG = x[3],
                                 ψ = x[4],
#                    α = x[7])* # Circular feed
                                 α = deg2rad(45.0))* # Circular feed
                      müller_sky(deg2rad(PAaz))*
                      lin_pol_source(flux = x[5],
                                     p = x[6],
                                     pa = x[7]))[i]
    for PAaz in PArange, i in 1:4]
    p1 = plot(PArange, pmodel[:,1], title = "I")
    scatter!(p1, pa, I, legend = :none, markersize = 3)
    p2 = plot(PArange, pmodel[:,2],title = "Q")
    scatter!(p2, pa, Q, legend = :none, markersize = 3)
    p3 = plot(PArange, pmodel[:,3],title = "U")
    scatter!(p3, pa, U, legend = :none, markersize = 3)
    p4 = plot(PArange, pmodel[:,4],title = "V")
    scatter!(p4, pa, V, legend = :none, markersize = 3)
    p = plot(p1, p2, p3, p4)
    display(p)
    @show sumsq, x
    return sumsq
end

# Set up NLopt
opt = Opt(:LN_COBYLA, 7)
ftol_rel!(opt, 0.000001)
#ftol_abs!(opt, 0.0004)

initial_step!(opt, [0.01, 0.01, 0.01, 0.01, 0.001, 0.05, 0.01])
#initial_step!(opt, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
lower_bounds!(opt, [-pi, -pi, -10, -pi, 0.01, 0.01, -2pi])
upper_bounds!(opt, [pi,   pi,  10,  pi, 0.05, 1.0,   2pi])
# initial_step!(opt, [0.1, 0.1, 0.1, 0.1, 0.1, 0.01])
# lower_bounds!(opt, [-pi, -pi, -10, -pi, 0.0, 0.01])
# upper_bounds!(opt, [pi, pi, 10, pi, 10.0, 1.0])

min_objective!(opt, (x, grad) -> cost_circ(x, grad)) #
#[-0.0024250372357195495,0.003151796261676011,-0.002372104102021719,0.11248945567718922,0.06121446041165875,0.01,-1.1219476970807971]
#minf, fit_param, ret = optimize!(opt, [0.0, 0.0, 0.0, 0.0, 0.1, 0.03, deg2rad(-64.0)])
minf, fit_param, ret = optimize!(opt, [-0.002, 0.003, -0.002, 0.11, 0.0225, 0.016, deg2rad(110.0)])
#minf, fit_param, ret = optimize!(opt, [0.0, 0.0, 0.0, 0.0, 0.1, 0.03])

model = Float64[(müller_rec(ϵ = fit_param[1],
                    ϕ = fit_param[2],
                    ΔG = fit_param[3],
                    ψ = fit_param[4],
#                    α = fit_param[7])* # Circular feed
                    α = deg2rad(45.0))* # Circular feed
         müller_sky(deg2rad(PAaz))*
         lin_pol_source(flux = fit_param[5],
                        p = fit_param[6],
                        pa = fit_param[7]))[i]
         for PAaz in PArange, i in 1:4]

# plot(plot(PArange, model[:,1], label = "I"), scatter!(pa, I, markersize = 3),
#      plot(PArange, model[:,2],label = "V"), scatter!(pa, V, markersize = 3),
#      plot(PArange, model[:,3],label = "U"), scatter!(pa, U, markersize = 3),
#      plot(PArange, -model[:,4],label = "Q"), scatter!(pa, Q, markersize = 3))

p1 = plot(PArange, model[:,1], title = "I")
scatter!(p1, pa, I, legend = :none, markersize = 3)
p2 = plot(PArange, model[:,2],title = "V")
scatter!(p2, pa, V, legend = :none, markersize = 3)
p3 = plot(PArange, model[:,3],title = "U")
scatter!(p3, pa, U, legend = :none, markersize = 3)
p4 = plot(PArange, -model[:,4],title = "Q")
scatter!(p4, pa, Q, legend = :none, markersize = 3)
plot(p1, p2, p3, p4)

max_stars_repo_name: KajWiik/PolarimetryCalibration

max_issues_repo_name: KajWiik/PolarimetryCalibration

max_forks_repo_name: KajWiik/PolarimetryCalibration


text:
function parse_string(str_in)
    arr_out = Array{Tuple{String, Int}, 1}() # initializ an array of tuples
    str_len = length(str_in)
    first_idx = 0 # keep track of the first index of the words
    tuple_idx = 1
    for i = 1:str_len # loop over every symbol in string
        if (first_idx < 1) # first_idx is reset to 0 after we've found a space
            first_idx = i
        end
        if (str_in[i] == ' ')
            push!(arr_out, (string(str_in[first_idx:i]),first_idx))
            first_idx = 0
            tuple_idx = tuple_idx + 1
            continue
        elseif i == str_len
            push!(arr_out, (string(str_in[first_idx:(i)]),first_idx))
        end
    end
    return arr_out
end

#struct Node
#    children::Dict{Char,Node}
#    posi::Array{Int}
#end

function build(list)
    toppNode = Node(Dict(),[])
    for (string, position) in list
        node = toppNode
        for char in string
            if !haskey(node.children, char)
                node.children[char] = Node(Dict(),[])
            end
            node = node.children[char]
        end
        push!(node.posi, position)
    end
    return toppNode
end

function positions(word, node, idx=1)
    if (idx > length(word))
        position = node.posi
    elseif (word[idx] == '?')
        position = []
        for key in collect(keys(node.children))
            append!(position, positions(word, node.children[key], idx + 1))
        end
    elseif (haskey(node.children, word[idx]))
        position = positions(word, node.children[word[idx]], idx + 1)
    else
        position = []
    end
    return position
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
using Base.Test
using DSGE

include(joinpath(dirname(@__FILE__()),"util.jl"))

my_tests = [
            "core",
            "parameters",
            "solve/gensys",
            "solve/solve",
            "estimate/kalman",
            "estimate/posterior",
            "estimate/hessizero",
            "estimate/hessian",
            "estimate/csminwel",
            "estimate/optimize",
            "estimate/eig",
            "estimate/metropolis_hastings",
            "models/m990/m990",
            "models/smets_wouters/smets_wouters",
            "data/misc",
            "data/load_data"
            ]

for test in my_tests
    test_file = string("$test.jl")
    @printf " * %s\n" test_file
    include(test_file)
end


text:
##############
# Constructors
##############

# Create some uncertain data of different types
o1 = UncertainValue(Normal, 0, 0.5)
o2 = UncertainValue(Normal, 2, 0.3)
o3 = UncertainValue(Uniform, 0, 4)
o4 = UncertainValue(Uniform, rand(100))
o5 = UncertainValue(rand(400))
o7 = CertainValue(2)
o8 = UncertainValue([2, 3, 4], [4, 5, 2])
o9 = UncertainValue([2, 4, 5, 2], rand(4))

uvals = [o1, o2, o3, o4, o5, o7, o8, o9]
UV = UncertainValueDataset(uvals)
UI = UncertainIndexDataset(uvals)
CUV = constrain(UV, TruncateQuantiles(0.1, 0.9))
CUI = constrain(UV, TruncateQuantiles(0.1, 0.9))

# From scalar vectors
@test UncertainIndexValueDataset(rand(30), rand(30)) isa UncertainIndexValueDataset

# Vectors
@test UncertainIndexValueDataset(uvals, uvals) isa UncertainIndexValueDataset

# Non-constrained datasets
@test UncertainIndexValueDataset(uvals, UV) isa UncertainIndexValueDataset
@test UncertainIndexValueDataset(UI, uvals) isa UncertainIndexValueDataset
@test UncertainIndexValueDataset(UI, UV) isa UncertainIndexValueDataset

# Constrained datasets
@test UncertainIndexValueDataset(uvals, CUV) isa UncertainIndexValueDataset
@test UncertainIndexValueDataset(CUI, uvals) isa UncertainIndexValueDataset
@test UncertainIndexValueDataset(CUI, CUV) isa UncertainIndexValueDataset

UIV = UncertainIndexValueDataset(UI, UV)

###########
# Iteration
###########
@test length(UIV) == 3
@test length([x for x in UIV]) == 3
@test UIV[1] isa Tuple{<:AbstractUncertainValue, <:AbstractUncertainValue}

###########
# Indexing
###########
@test UIV[1] isa Tuple{<:AbstractUncertainValue, <:AbstractUncertainValue}
@test UIV[end] isa Tuple{<:AbstractUncertainValue, <:AbstractUncertainValue}
@test UIV[1:end] isa AbstractVector
@test UIV[[1, 2]] isa AbstractVector
@test UIV[:] isa AbstractVector

@test index(UIV, 1) isa AbstractUncertainValue
@test index(UIV, 1:2) isa AbstractVector{<:AbstractUncertainValue}