text:
using Base.Test
using Distributions
using Klara

ρ = Hyperparameter(:ρ, 1)

p1 = BasicContUnvParameter(
  :p1,
  2,
  signature=:low,
  setpdf=(state, states) -> Normal(states[1].value*states[3].value, sqrt(1-abs2(states[1].value)))
)

p2 = BasicContUnvParameter(
  :p2,
  3,
  signature=:low,
  setpdf=(state, states) -> Normal(states[1].value*states[2].value, sqrt(1-abs2(states[1].value)))
)

model = GenericModel([ρ, p1, p2], [ρ p1; ρ p2; p1 p2])

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

vstate = [BasicUnvVariableState(0.8), BasicContUnvParameterState(5.1), BasicContUnvParameterState(2.3)]

outopts = [Dict(:monitor=>[:value]), Dict(:monitor=>[:value])]

job = BasicGibbsJob(
  model,
  [2, 3],
  [nothing, nothing],
  mcrange,
  vstate,
  outopts,
  true,
  false,
  false
)

@time run(job)

# output(job)

chains = Dict(job)

[mean(chains[k]) for k in keys(chains)]

cor(chains[:p1].value, chains[:p2].value)


text:
function __init__()
    global UA_TYPES
    global UA_TYPES_PTRS
    global UA_TYPES_MAP
    UA_TYPES[] = cglobal((:UA_TYPES, libopen62541), UA_DataType)

    for i = 0:UA_TYPES_COUNT-1
        UA_TYPES_PTRS[i] = UA_TYPES[] + sizeof(UA_DataType) * i

        typename = "UA_" * unsafe_string(unsafe_load(UA_TYPES_PTRS[i]).typeName)
        UA_TYPES_MAP[i+1] = eval(Symbol(typename))
    end
end

ext: jl


text:
## header to provide surface and context
using Cairo
using Printf

c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

set_font_face(cr, "Sans 16")

text(cr,40.0,40.0,"EF",markup=true)
text(cr,40.0,60.0,"EGDC")

text(cr,40.0,90.0,"EFGH",markup=true)
text(cr,40.0,110.0,"EGDC",markup=true)

text(cr,40.0,140.0,"GEFJ",markup=true)
text(cr,40.0,160.0,"EGDC")

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_pango_text1.png");


text:
@testset "context" begin

@test ctx == CuCurrentContext()
@test ctx === CuCurrentContext()

let ctx2 = CuContext(dev)
    @test ctx2 == CuCurrentContext()    # ctor implicitly pushes
    activate(ctx)
    @test ctx == CuCurrentContext()

    @test device(ctx2) == dev

    CUDAdrv.unsafe_destroy!(ctx2)
end

let global_ctx2 = nothing
    CuContext(dev) do ctx2
        @test ctx2 == CuCurrentContext()
        @test ctx != ctx2
        global_ctx2 = ctx2
    end
    @test !CUDAdrv.isvalid(global_ctx2)
    @test ctx == CuCurrentContext()

    @test device(ctx) == dev
    @test device() == dev
    synchronize()
end

end

@testset "primary context" begin

pctx = CuPrimaryContext(dev)

@test !isactive(pctx)
unsafe_reset!(pctx)
@test !isactive(pctx)

@test flags(pctx) == CUDAdrv.CTX_SCHED_AUTO
setflags!(pctx, CUDAdrv.CTX_SCHED_BLOCKING_SYNC)
@test flags(pctx) == CUDAdrv.CTX_SCHED_BLOCKING_SYNC

let global_ctx = nothing
    CuContext(pctx) do ctx
        @test CUDAdrv.isvalid(ctx)
        @test isactive(pctx)
        global_ctx = ctx
    end
    @test !isactive(pctx)
    @test !CUDAdrv.isvalid(global_ctx)
end

CuContext(pctx) do ctx
    @test CUDAdrv.isvalid(ctx)
    @test isactive(pctx)

    unsafe_reset!(pctx)

    @test !isactive(pctx)
    @test !CUDAdrv.isvalid(ctx)
end

let
    @test !isactive(pctx)

    ctx1 = CuContext(pctx)
    @test isactive(pctx)
    @test CUDAdrv.isvalid(ctx1)

    unsafe_reset!(pctx)
    @test !isactive(pctx)
    @test !CUDAdrv.isvalid(ctx1)
    CUDAdrv.valid_contexts

    ctx2 = CuContext(pctx)
    @test isactive(pctx)
    @test !CUDAdrv.isvalid(ctx1)
    @test CUDAdrv.isvalid(ctx2)

    unsafe_reset!(pctx)
end

end

@testset "cache config" begin

cache_config!(CUDAdrv.FUNC_CACHE_PREFER_L1)
@test cache_config() == CUDAdrv.FUNC_CACHE_PREFER_L1

end

@testset "shmem config" begin

shmem_config!(CUDAdrv.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE)
@test shmem_config() == CUDAdrv.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE

end

@testset "limits" begin

lim = limit(CUDAdrv.LIMIT_DEV_RUNTIME_SYNC_DEPTH)

lim += 1
limit!(CUDAdrv.LIMIT_DEV_RUNTIME_SYNC_DEPTH, lim)
@test lim == limit(CUDAdrv.LIMIT_DEV_RUNTIME_SYNC_DEPTH)

end


text:
include("IOM/src/share/LogSystem.jl")
include("IOM/src/share/PolelikeCoordinate.jl")

include("IOM/src/models/EMOM/ENGINE_EMOM.jl")
include("IOM/src/driver/driver_working.jl")

using MPI
using CFTime
using Dates
using ArgParse
using .PolelikeCoordinate
using .LogSystem

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin
        "--read-restart"
            help = "Configuration file."
            arg_type = Bool

        "--stop-n"
            help = "Core of the model."
            arg_type = Int64

        "--time-unit"
            help = "."
            arg_type = String
    end

    return parse_args(s)
end

parsed = parse_commandline()

t_start = DateTimeNoLeap(1, 1, 1)

time_unit = Dict(
    "year" => Dates.Year,
    "month" => Dates.Month,
    "day"   => Dates.Day,
)[parsed["time-unit"]]

t_simulation = time_unit(parsed["stop-n"])
t_end = t_start + t_simulation
MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
Δt = Second(86400)

include("config.jl")

if rank == 0
    gf = PolelikeCoordinate.CurvilinearSphericalGridFile(
        config["MODEL_CORE"]["domain_file"],
        R  = 6371229.0,
        Ω  = 2π / (86400 / (1 + 365/365)),
    )

end

coupler_funcs = (
    master_before_model_init = function()
        return parsed["read-restart"], t_start
    end,

    master_after_model_init! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] After model init")

    end,
    master_before_model_run! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] Before model run")
        writeLog("[Coupler] This is where flux exchange happens.")

        global comm, rank, t_end

        is_master = rank == 0

        # Test if t_end is reached
        return_values = nothing
        if is_master
            t_end_reached = OMDATA.clock.time >= t_end

            if t_end_reached
                return_values = ( :END, Δt, t_end_reached )
            else
                return_values = ( :RUN, Δt, t_end_reached )
            end
        end

        # Deal with coupling
        if is_master

            writeLog("[Coupler] Need to broadcast forcing fields.")
            # compute flux
            lat = gf.yc
            lon = gf.xc

#            OMDATA.x2o["SWFLX"][1, :, :] .= - (cos.(deg2rad.(lat).+0.1) .+ 1) / 2 .* (sin.(deg2rad.(lon)) .+ 1)/2 * 100.0
#            OMDATA.x2o["TAUX_east"][1, :, :]   .= 0.2 * (cos.(deg2rad.(lat).+0.1) .+ 1) / 2
#            OMDATA.x2o["TAUY_north"][1, :, :]  .= 0.1 * (sin.(deg2rad.(lon)) .+ 1) / 2
#cos.(deg2rad.(lat))

        end

        return return_values

    end,

    master_after_model_run! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] After model run")
    end,

    master_finalize! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] Finalize")

    end
)

runModel(
    ENGINE_EMOM,
    coupler_funcs,
    config,
)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Chemostat_Rath2017
using Test

@testset "Chemostat_Rath2017.jl" begin
    # Write your own tests here.
    @test true
end


text:
# Autogenerated wrapper script for TMatrix_jll for x86_64-w64-mingw32-libgfortran5
export tmatrix_fixed_orientation, tmatrix_random_orientation

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("TMatrix")
JLLWrappers.@declare_library_product(tmatrix_fixed_orientation, "libtmatrixfo.dll")
JLLWrappers.@declare_library_product(tmatrix_random_orientation, "libtmatrixro.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        tmatrix_fixed_orientation,
        "bin\\libtmatrixfo.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        tmatrix_random_orientation,
        "bin\\libtmatrixro.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# -------------------------------------------
# LocalProjections:
# For each individual tide-gauge location:
# - Read local observations
# - Compute local extrapolation
# - Read local NCA5 projections
# - Put everything on a year 2000 baseline
# Save all to NOAA_TR_local_projections.nc
# -------------------------------------------

module LocalProjections

using NetCDF
using DelimitedFiles
using Plots
using Statistics
using NCDatasets
using MAT
using Interpolations
using CSV

dir_code = homedir()*"/Projects/2021_NOAA_TR/Code/"
include(dir_code*"Masks.jl")
include(dir_code*"Hector.jl")
include(dir_code*"RegionalProjections.jl")
include(dir_code*"ProcessObservations.jl")

function RunLocalProjections(settings)
    println("\nLocal observations, trajectories and scenarios...")
    local_obs  = ProcessObservations.ReadLocalObs(settings)
    ExtrapolateLocalTrajectory!(local_obs,settings)
    NCA5_local,years_NCA5,pct_NCA5 = ReadLocalProjections(local_obs,settings)
    correct_baselines!(local_obs,NCA5_local, settings)
    save_data(local_obs,NCA5_local, years_NCA5, pct_NCA5,settings)
    println("Local observations, trajectories and scenarios done\n")
end

function ExtrapolateLocalTrajectory!(local_obs,settings)
    # Compute and extrapolate the trajectory in observed sea level
    println("  Extrapolating tide-gauge records to compute trajectory...")
    # Prepare design matrices
    y_acc = findall(in(settings["years_trajectory"]),settings["years_tg"])
    amat = ones(length(y_acc),3)
    amat[:,2] = settings["years_tg"][y_acc] .- mean(settings["years_tg"][y_acc])
    amat[:,3] = (settings["years_tg"][y_acc] .- mean(settings["years_tg"][y_acc])).^2
    amat_extend = ones(length(settings["years_trajectory"]),3)
    amat_extend[:,2] .= settings["years_trajectory"] .- mean(settings["years_tg"][y_acc])
    amat_extend[:,3] = (settings["years_trajectory"] .- mean(settings["years_tg"][y_acc])).^2

    # Do the extrapolation
    local_obs["obs_length"] = zeros(Int32,length(local_obs["name"])); # Number of observation-years used to estimate trajectory
    local_obs["obs_lt_30"] = zeros(Int32,length(local_obs["name"])); # Number of observation-years used to estimate trajectory
    local_obs["rsl_trend"] = zeros(Float32,length(local_obs["name"]),3);
    local_obs["rsl_accel"] = zeros(Float32,length(local_obs["name"]),3);
    local_obs["rsl_trajectory"] = zeros(Float32,length(local_obs["name"]),length(settings["years"]),3);
    for tg in 1:length(local_obs["name"])
        local_obs["obs_length"][tg] = sum(isfinite.(local_obs["rsl"][tg,y_acc]))
        local_obs["obs_length"][tg] >= 30 ? local_obs["obs_lt_30"][tg] = 1 : local_obs["obs_lt_30"][tg] = 0
        trend_file = Hector.EstTrend(settings["years_tg"][y_acc],local_obs["rsl"][tg,y_acc];accel=true,model="Powerlaw",SA=false,SSA=false,monthly=false,tref=mean(settings["years_tg"][y_acc]))
        μ_sol = [trend_file["bias"],trend_file["trend"],trend_file["accel"]]
        σ_sol= [trend_file["bias_sigma"],trend_file["trend_sigma"],trend_file["accel_sigma"]]
        sol_arr = reshape(μ_sol,(1,3)) .+ randn(5000,3) .* reshape(σ_sol,(1,3))
        trajectory_arr = zeros(Float32,length(settings["years_trajectory"]),5000)
        [trajectory_arr[:,i] = amat_extend * @views sol_arr[i,:] for i ∈ 1:5000]
        trajectory_stats = zeros(Float32,length(settings["years_trajectory"]),3)
        [trajectory_stats[t,:] = quantile((@views trajectory_arr[t,:]),[0.17,0.50,0.83]) for t ∈ 1:length(settings["years_trajectory"])]
        local_obs["rsl_trajectory"][tg,:,:] = LinearInterpolation((settings["years_trajectory"],[1.0f0:3.0f0...]),trajectory_stats,extrapolation_bc=NaN32)[settings["years"],[1.0f0:3.0f0...]]
        local_obs["rsl_trend"][tg,:] = [μ_sol[2]-σ_sol[2],μ_sol[2],μ_sol[2]+σ_sol[2]]
        local_obs["rsl_accel"][tg,:] = 2 .*[μ_sol[3]-σ_sol[3],μ_sol[3],μ_sol[3]+σ_sol[3]]
    end
    return nothing
end

function ReadLocalProjections(local_obs,settings)
    # Locations and time
    fn_grid = settings["dir_NCA5"]*"NCA5_Low_grid.nc"
    lon_NCA5 = ncread(fn_grid,"lon")
    lat_NCA5 = ncread(fn_grid,"lat")
    years_NCA5 = convert.(Float32,ncread(fn_grid,"years",start=[1],count=[-1]))
    pct_NCA5 = convert.(Float32,ncread(fn_grid,"percentiles"))

    fh_tg = Dataset(settings["dir_NCA5"]*"NCA5_Low_tg.nc","r")
    psmsl_loc = fh_tg["loc"][:]
    close(fh_tg)

    NCA5_loc_is_present = zeros(Bool,length(local_obs["name"]))
    NCA5_loc_tg = zeros(Int32,length(local_obs["name"]))
    NCA5_loc_grid = zeros(Int,length(local_obs["name"]),2)
    for tg in 1:length(local_obs["name"])
        loc_idx = findfirst(==(local_obs["psmsl_id"][tg]),psmsl_loc)
        if ~isnothing(loc_idx)
            # Tide gauge is in PSMSL/local database: use tide-gauge projection
            NCA5_loc_is_present[tg] = true
            NCA5_loc_tg[tg] = loc_idx
            NCA5_loc_grid[tg,1] = argmin(@. abs(local_obs["coords"][tg,1]-lon_NCA5))
            NCA5_loc_grid[tg,2] = argmin(@. abs(local_obs["coords"][tg,2]-lat_NCA5))
        else
            # Tide gauge is not in PSMSL/local database: sample location from the grid
            NCA5_loc_is_present[tg] = false
            NCA5_loc_tg[tg] = -1
            NCA5_loc_grid[tg,1] = argmin(@. abs(local_obs["coords"][tg,1]-lon_NCA5))
            NCA5_loc_grid[tg,2] = argmin(@. abs(local_obs["coords"][tg,2]-lat_NCA5))
        end
    end

    # Create data sctructure
    NCA5_local = Array{Dict}(undef,length(local_obs["name"]))
    for tg in 1:length(local_obs["name"])
        NCA5_local[tg] = Dict()
        NCA5_local[tg]["Loc_is_in_proj"] = NCA5_loc_is_present[tg]
        for scn in settings["NCA5_scenarios"]
            NCA5_local[tg][scn] = Dict()
        end
    end

    # Read data
    for scenario in settings["NCA5_scenarios"]
        println("   Scenario "*scenario*"...")
        fn_grid = settings["dir_NCA5"]*"NCA5_"*scenario*"_grid.nc"
        fn_tg = settings["dir_NCA5"]*"NCA5_"*scenario*"_tg.nc"
        for prc in settings["processes"]
            NCA5_prc_grid = convert.(Float32,ncread(fn_grid,prc));
            NCA5_prc_tg = convert.(Float32,ncread(fn_tg,prc));
            for tg in 1:length(local_obs["name"])
                if NCA5_loc_is_present[tg] # Use TG scenario
                    NCA5_local[tg][scenario][prc] = NCA5_prc_tg[NCA5_loc_tg[tg],:,:]
                else # Use gridded scenario
                    NCA5_local[tg][scenario][prc] = NCA5_prc_grid[NCA5_loc_grid[tg,1],NCA5_loc_grid[tg,2],:,:]
                end
            end
        end
    end
    return NCA5_local,years_NCA5,pct_NCA5
end

function correct_baselines!(local_obs,NCA5_local, settings)
    println("  Correcting baselines for trajectory, observations and scenarios...")
    # Correct baseline: baseline for everything is trajectory value in 2000
    traj_baseline = findall(in(settings["years_baseline"]),settings["years"])
    traj_nca5_start = findfirst(settings["years"].==2005)

    for tg in 1:length(local_obs["name"])
        # 1. Baseline trajectory
        local_obs["rsl_trajectory"][tg,:,:] .-= local_obs["rsl_trajectory"][tg,traj_baseline,2]

        # 2. Baseline observations
        acc_obs = isfinite.(local_obs["rsl"][tg,:])
        yr_int = intersect(settings["years_tg"][acc_obs],settings["years_trajectory"])
        obs_mn = mean(local_obs["rsl"][tg,(findall(in(yr_int),settings["years_tg"]))])
        traj_mn = mean(local_obs["rsl_trajectory"][tg,(findall(in(yr_int),settings["years"])),2])
        local_obs["rsl"][tg,:] = local_obs["rsl"][tg,:] .- obs_mn .+ traj_mn

        # 3. Baseline scenarios
        for scenario in settings["NCA5_scenarios"]
            NCA5_local[tg][scenario]["total"] .+= local_obs["rsl_trajectory"][tg,traj_nca5_start,2]
            NCA5_local[tg][scenario]["verticallandmotion"] .+= local_obs["rsl_trajectory"][tg,traj_nca5_start,2]
        end
    end
    return nothing
end

function save_data(local_obs,NCA5_local, years_NCA5, pct_NCA5,settings)
    println("  Saving data...");
    scn_proj = Dict()
    for scenario ∈ settings["NCA5_scenarios"]
        scn_proj[scenario] = Dict()
        for prc in settings["processes"]
            scn_proj[scenario][prc] = Array{Float32}(undef,length(local_obs["name"]),length(settings["years"]),3)
            for tg in 1:length(local_obs["name"])
                scn_proj[scenario][prc][tg,:,:] = LinearInterpolation((years_NCA5, pct_NCA5), NCA5_local[tg][scenario][prc],extrapolation_bc=NaN32)(settings["years"],pct_NCA5)
            end
        end
    end

    rsl_obs = zeros(Float32,length(local_obs["name"]),length(settings["years"])) .* NaN32
    acc_obs = findall(in(settings["years_tg"]),settings["years"])
    rsl_obs[:,acc_obs] = local_obs["rsl"]

    fh = Dataset(settings["fn_proj_lcl"],"c")
    fh.attrib["title"] = "Local projections"
    fh.attrib["description"] = "Local projections, trajectory, and observations at tide-gauge locations for the interagency report: Global and Regional Sea Level Rise Scenarios for the United States: Updated Mean Projections and Extreme Water Level Probabilities Along U.S. Coastlines"
    fh.attrib["processes"] = "AIS: Antarctic Ice Sheet, GIS: Greenland, glaciers: Glaciers and Ice Caps, landwaterstorage: Liquid water storage changes on land, oceandynamics: Ocean dynamics and global thermosteric expansion, verticallandmotion: Vertical land motion, total: All processes combined."

    defDim(fh,"years", length(settings["years"]))
    defDim(fh,"percentiles",3)
    defDim(fh,"tg", length(local_obs["name"]))

    defVar(fh,"tg",local_obs["name"],("tg",),deflatelevel=5, attrib = Dict("description" => "Name of the tide gauge", "units" => "-"))
    defVar(fh,"lon",local_obs["coords"][:,1],("tg",),deflatelevel=5, attrib = Dict("description" => "Longitude of the tide gauge", "units" => "degrees East"))
    defVar(fh,"lat",local_obs["coords"][:,2],("tg",),deflatelevel=5, attrib = Dict("description" => "Latitude of the tide gauge", "units" => "degrees North"))
    defVar(fh,"years",settings["years"],("years",),deflatelevel=5, attrib = Dict("units" => "years"))
    defVar(fh,"percentiles",convert.(Int32,[17,50,83]),("percentiles",),deflatelevel=5, attrib = Dict("description" => "The percentile of the projection. The 50th percentile shows the median projection, and the 17th and 83rd show the upper- and lower bound on the 1 sigma level.", "units" => "-"))

    defVar(fh,"PSMSL_id",convert.(Int32,local_obs["psmsl_id"]),("tg",),deflatelevel=5, attrib = Dict("description" => "PSMSL ID to which the tide gauge corresponds. Stations that are not in PSMSL are flagged -1."))
    defVar(fh,"QC_flag",convert.(Int8,local_obs["qc_flag"]),("tg",),deflatelevel=5, attrib = Dict("description" => "Quality control flag. Value of 1 means possible issue. Currently manual process."))
    # Write trajectory
    defVar(fh,"rsl_obs",rsl_obs,("tg","years"),deflatelevel=5, attrib = Dict("description" => "Observed local sea-level changes, relative to year 2000 (tide-gauge observations).", "units" => "mm"))
    defVar(fh,"number_of_observation_years_for_trajectory_estimation",local_obs["obs_length"],("tg",),deflatelevel=5, attrib = Dict("description" => "Number of observation-years used to estimate the trajectory.", "units" => "-"))
    defVar(fh,"number_of_observation_years_for_trajectory_estimation_geq_30",local_obs["obs_lt_30"],("tg",),deflatelevel=5, attrib = Dict("description" => "Flag if number of observation-years used for the trajectory exceeds 30. 1: 30 years or more, 0: less than 30 years.", "units" => "-"))
    defVar(fh,"rsl_trajectory",local_obs["rsl_trajectory"],("tg","years","percentiles"),deflatelevel=5, attrib = Dict("description" => "Estimated local sea-level trajectory, relative to year 2000.", "units" => "mm"))
    defVar(fh,"rsl_trend",local_obs["rsl_trend"],("tg","percentiles"),deflatelevel=5, attrib = Dict("description" => "Linear trend of estimated trajectory", "units" => "mm yr-1"))
    defVar(fh,"rsl_accel",local_obs["rsl_accel"],("tg","percentiles"),deflatelevel=5, attrib = Dict("description" => "Acceleration (half quadratic) of estimated trajectory", "units" => "mm yr-2"))
    # Write scenarios
    for scn ∈ settings["NCA5_scenarios"]
        for prc in settings["processes"]
            defVar(fh,"rsl_"*prc*"_"*scn,scn_proj[scn][prc],("tg","years","percentiles",),deflatelevel=5, attrib = Dict("description" => "Projected local sea level for scenario "*scn*" and process "*prc*", relative to year 2000.", "units" => "mm"))
        end
    end
    close(fh)
    return nothing
end

function example_plots(local_obs,NCA5_local, years_NCA5, pct_NCA5,settings)
    # Make a simple plot of all local time series
    for tg in 1:length(local_obs["name"])
        scn_cl = cgrad(:viridis, 3, categorical = true)
        plot(settings["years_tg"],local_obs["rsl"][tg,:],color=:red,linewidth=2)
        plot!(settings["years"],local_obs["rsl_trajectory"][tg,:,2],color=:blue,linewidth=2)
        plot!(years_NCA5,NCA5_local[tg]["Low"]["total"][:,2],color=scn_cl[1],linewidth=2)
        plot!(years_NCA5,NCA5_local[tg]["Int"]["total"][:,2],color=scn_cl[2],linewidth=2)
        plot!(years_NCA5,NCA5_local[tg]["High"]["total"][:,2],color=scn_cl[3],linewidth=2,legend=false)
        xlims!((1970,2050))
        ylims!((-200,1000))
        title!(local_obs["name"][tg])
        savefig(homedir()*"/Scratch/Plots/"*string(tg)*".png")
    end
end

end


text:
println(-42)
println(40+2)
println(44-2)
println(6*7)
println(div(84,9))
println(div(-84,9))
println(div(84,-9))
println(div(-84,-9))
println(84 % 9)
println(-84 % 9)
println(84 % -9)
println(-84 % -9)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
const TYPES = (Bool, Int32, Int64, Int128, Float32, Float64, String, UInt8)

read_column(path, col_num) = read_column(path, metadata(path), col_num)

function read_column(path, filemetadata, col_num)
    par = ParFile(path)

    T = TYPES[filemetadata.schema[col_num+1]._type+1]
    # TODO detect if missing is necessary
    res = Vector{Union{Missing, T}}(missing, nrows(par))
    write_cursor = 1
    for row_group in filemetadata.row_groups
        pgs = pages(par, row_group.columns[col_num])

        drop_page_count = 0
        # is the first page a dictionary page
        # this is not the case for boolean values for example
        if isfilled(pgs[1].hdr, :dictionary_page_header)
            # the first page is almost always the dictionary page
            dictionary_page = pgs[1]
            drop_page_count = 1
            dictionary_of_values = T.(values(par, dictionary_page)[1])
        end

        # TODO deal with other types of pages e.g. dataheaderv2

        # everything after the first data datapages
        for data_page in Base.Iterators.drop(pgs, drop_page_count)
            vals, definitions, decode = values(par, data_page)

            @assert all(in((0, 1)), definitions)

            l = sum(==(1), definitions)
            # if all definitions values are 1 then it's not used
            definitions_not_used = all(==(1), definitions)

            # data_page can be either
            # * dictionary-encoded in which case we should look into the dictionary
            # * plained-encoded in which case just return the values
            page_encoding = Diban.page_encoding(data_page)

            if page_encoding == Encoding.PLAIN_DICTIONARY
                if definitions_not_used
                    res[write_cursor:write_cursor+l-1] .= dictionary_of_values[vals.+1]
                else
                    val_index = 1
                    for (offset, definitions) in enumerate(definitions)
                        if definition != 0
                            value = vals[val_index]
                            res[write_cursor+offset-1] = dictionary_of_values[value + 1]
                            val_index += 1
                        end
                    end
                end
            elseif page_encoding == Encoding.PLAIN
                if definitions_not_used
                    res[write_cursor:write_cursor+l-1] .= T.(vals)
                else
                    val_index = 1
                    for (offset, definition)  in enumerate(definitions)
                        if definition != 0
                            value = vals[val_index]
                            res[write_cursor+offset-1] = T(value)
                            val_index += 1
                        end
                    end
                end
            else
                error("page encoding not supported yet")
            end

            write_cursor += length(definitions)
        end
    end
    return res
end

()


text:
module tele
using NCDatasets, Dates, PyCall, Revise, PyPlot, Statistics, NaNMath

export read_netcdf, contplot, anomaly, batch_open, character,filtseriesplot,
normalequation, fftplot, yearly_mean, spatial_mean, temp_mean, plot_eofs, plot_corr, plot_reg

const ccrs = PyNULL()
const cfeature = PyNULL()
const mpatches = PyNULL()
const ticker = PyNULL()

"""
   function nc_dims(nc)
   Get all important netcdf dims assume CF-conventions
# Arguments
- `nc`: NCDataset
# Output
- `lat`: all lat values
- `lon`: all lon values
- `time`: time values (already converted to Datetime)
"""
function nc_dims(nc)
   lat = getindex(nc, "lat")
   lon = getindex(nc, "lon")
   time = getindex(nc, "time")
   if length(size(lat)) > 1
      a = lat[1,:]
      o = lon[:,1]
   else
      a = lat
      o = lon
   end

   return a, o, time
end

"""
   function read_netcdf(path, varname, bbox)
   Read in NetCDF file and spatially subset into square bounding box
# Arguments#
- `path`: string, path to netcdf file
- `varname`: string, name of variable you want to read
- `bbox`: geographic bounding box in format [latmin, latmax, lonmin, lonmax]
# Outputlat_diff = pi*r_e / 180 #km/°latlat_diflat_diff = pi*r_e / 180 #km/°lat
f = pi*r_e / 180 #km/°lat

- `lat`: subsetted lat
- `lon`: subsetted lon
- `time`: subsetted time
- `data`: subsetted data
"""
function read_netcdf(path, varname, bbox)
   nc = NCDataset(path, "r")
   lat, lon, time = nc_dims(nc)
   #find lat lon box indices
   lat_flag = lat[begin] > lat[end] ? true : false #check if lat is reversed
   a1 = searchsortedfirst(lat, bbox[1], rev = lat_flag)
   a2 = searchsortedfirst(lat, bbox[2], rev = lat_flag)
   o1 = searchsortedfirst(lon, bbox[3])
   o2 = searchsortedfirst(lon, bbox[4])
   println([a1,a2,o1,o2])

   #is our longitude array continuous - not true for CESM
   if maximum(abs.(diff(lon))) > 100
      ind = findall(>(100), abs.(diff(lon)))[1]
      lon1 = lon[begin:ind-1]
      lon2 = lon[ind:end]

      o1 = searchsortedfirst(lon1, bbox[3])

      bool = zeros(2)
      if o1 == 1
         o1 = searchsortedfirst(lon2, bbox[3])
         bool[1] = 1
      end

      o2 = searchsortedfirst(lon1, bbox[4])
      if o2 == 1
         o2 = searchsortedfirst(lon2, bbox[4])
         bool[2] = 1
      end

      println("bool is " *string(bool))
      if bool == zeros(2)
         println("first bool : "* string([o1,o2]))
         lon_sub = lon1[minimum([o1,o2]):maximum([o1,o2])]
         subdata = nc[varname][minimum([o1,o2]):maximum([o1,o2]), minimum([a1, a2]):maximum([a1, a2]), :]

      elseif bool == [1,1]
         lon_sub = lon1[minimum([o1,o2]):maximum([o1,o2])]
         o1 = length(lon1) + o1
         o2 = length(lon1) + o2
         subdata = nc[varname][minimum([o1,o2]):maximum([o1,o2]), minimum([a1, a2]):maximum([a1, a2]), :]
      #FILLOWING CODE IS NOT TESTED - this is where I'd put my good code (if I had it)
      # elseif bool == [0, 1]
      #    lon_sub = vcat(lon1[o1:end], lon2[begin:o2])
      #    subdata1 = nc[varname][o1:length(lon1), minimum([a1, a2]):maximum([a1, a2]), :]
      #    subdata2 = nc[varname][length(lon1):o2+length(lon1), minimum([a1, a2]):maximum([a1, a2]), :]
      #    subdata = cat(subdata1, subdata2, dims = 2)

      #this works! (amazing)
      elseif bool == [1, 0]
         lon_sub = vcat(lon2[o1:end], lon1[begin:o2])
         subdata1 = nc[varname][length(lon1) + o1:end, minimum([a1, a2]):maximum([a1, a2]), :]
         subdata2 = nc[varname][begin:o2, minimum([a1, a2]):maximum([a1, a2]), :]
         subdata = cat(subdata1, subdata2, dims = 1)
      end

   else
      subdata = nc[varname][minimum([o1,o2]):maximum([o1,o2]), minimum([a1, a2]):maximum([a1, a2]), :]
      lon_sub = lon[minimum([o1,o2]):maximum([o1,o2])]
   end

   na = nc[varname].attrib["missing_value"]
   replace!(subdata,na => NaN)
   replace!(subdata,missing=> NaN)
   subdata = Float64.(subdata)
   return lat[minimum([a1, a2]):maximum([a1, a2])], lon_sub, time, subdata

end

#this is hacky!
function load_python()
   copy!(ccrs, pyimport("cartopy.crs"))
   copy!(cfeature, pyimport("cartopy.feature"))
   copy!(mpatches, pyimport("matplotlib.patches"))
   copy!(ticker, pyimport("matplotlib.ticker"))
end

"""
   function contplot(lat, lon, data, s)
   Make a contour plot of data (assumes SouthPolarStereo view) and save
# Arguments#
- `lat`
- `lon`
- `data`: data in (lat x lon x time) format
- `s`: file path to save to
"""
function contplot(lat,lon, data,s)
   load_python()
   figure()
   ax = subplot(projection = ccrs.PlateCarree())
   ax.gridlines()
   #ax.set_extent([-180, 0, 50, 90], ccrs.PlateCarree())
   gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=true,
                     linewidth=2, color="gray", alpha=0.5, linestyle="--")
   gl.right_labels = false
   ax.coastlines()
   max = NaNMath.maximum(vcat(data...))
   min = NaNMath.minimum(vcat(data...))
   println(min, max)
   val = maximum([abs(max), abs(min)])

   cf = contourf(lon,lat,transpose(data), cmap = "coolwarm", levels = range(-val, stop = val, length = 10), transform=ccrs.PlateCarree())
   cb = colorbar(cf, fraction = 0.015)
   cb.set_label("SST [σ]")
   tight_layout()
   savefig(s,bbox_inches="tight")
   close()
end

"""
   function anomaly(data, time)
   Remove seasonality by computing the anomaly (stdev away from mean) w.r.t. behavior in month
# Arguments#
- `data`: array of data in (lat x lon x time) format
- `time`: array of DateTime
# Output
- `anomaly`: array of same dimensions of `data` where each value is the anomaly in stdev
"""
function anomaly(data, time)
   anomaly = similar(data)

   #make list of months
   m = []
   for i in time
      push!(m, Dates.month(i))
   end

   #compute average for each month
   marray = Array{Float32}(undef, size(data)[1], size(data)[2],12)
   sarray = Array{Float32}(undef, size(data)[1], size(data)[2],12)
   for i in range(1, stop = 12, step = 1) #iterate over each month, collecting all vals

      #find indices of months of interest
      bool_t_ind = findfirst.(isequal.(m), (i,)).== nothing
      t_ind = findall(x->x==0, bool_t_ind)
      marray[:, :, i] .= temp_mean(data[:,:,t_ind])
      sarray[:,:, i] .= temp_std(data[:,:,t_ind])
   end

   for i in range(1, stop = length(time), step = 1)
      m_index = Dates.month(time[i])
      anomaly[:,:,i] = (data[:,:,i] .- marray[:,:, m_index]) ./ sarray[:,:,m_index]
   end

   return anomaly

end

"""
   function yearly_mean(time, vals)
   Compute annual average of 1d array (assumes data is ordered AND continuous)
# Arguments#
- `time`: array of time values (DateTime)
- `vals`: array of associated values
# Output
- `ann_avg`: annual average of `vals`
- `years`: list of years associated with `ann_avg`
"""
function yearly_mean(time, vals)
   initial_month = Dates.month(time[1])
   initial_year = Dates.year(time[1])
   stop_month = Dates.month(time[end])
   stop_year = Dates.year(time[end])

   years = range(initial_year, stop = stop_year, step = 1) #array of all possible years

   ann_avg = zeros(length(years))
   #compute first and last values, because we might not have 12 months in first or last year
   ann_avg[1] = mean(vals[1:13-initial_month])
   ann_avg[end] = mean(vals[end - stop_month+1:end])

   #iterate through, computing the average for every 12 months, starting from end of year 1
   for i in range(2, length(ann_avg)-1, step = 1)
      start_index = 12-initial_month + (i-1)*12
      stop_index = start_index + 12
      ann_avg[i] = mean(vals[start_index:stop_index])
   end

   return ann_avg, years
end

"""
   function character(series)
   Compute statistics for a series
# Arguments#
- `series`: array of floats, can have NaN values (will be ignored)
# Output
- `dict`: dictionary of statistics
"""
function character(series)
    min = NaNMath.minimum(series)
    max = NaNMath.maximum(series)
    mean = NaNMath.mean(series)
    median = NaNMath.median(series)
    std = NaNMath.std(series)
    dict = Dict(:min=>min, :max=> max, :mean=>mean, :median=>median, :std=>std)
    return dict
end

"""
   function filtseriesplot(series)
      This is a pretty specific function - don't need to go into everything
      Makes plot with 2 subplots, each has raw data, filtered data, and a trendline
"""
function filtseriesplot(series1, series1filt, series1time, series2, series2filt, series2time, t, s)
   figure(figsize = (10, 4))
   subplot(1,2,1)
   shift = 50
   plot(series1time, series1)
   plot(series1time[shift:end], series1filt[shift:end])
   x_ind1 = range(1, stop = length(series1time[shift:end]), step = 1)
   s1_avg = mean(series1filt[shift:end])
   series1filt_slope = normalequation(x_ind1, series1filt[shift:end] .- s1_avg)
   series1filt_ne = series1filt_slope*x_ind1 .+ s1_avg
   plot(series1time[shift:end], series1filt_ne, color = "r")
   xpos = series1time[1000]
   ypos = s1_avg
   legend(["CESM", "CESM smoothed", "CESM smoothed trend: "* string(round(series1filt_slope, digits =6))])
   xlabel("time")

   if s == "images/anom.pdf"
      ylabel("SST Anomaly [σ]")
   else
      ylabel("SST [°C]")
   end

   subplot(1,2,2)
   plot(series2time, series2)
   shift = 50
   plot(series2time[shift:end], series2filt[shift:end])
   x_ind2 = range(1, stop = length(series2time[shift:end]), step = 1)
   s2_avg = mean(series2filt[shift:end])
   print(s2_avg)
   series2filt_slope = normalequation(x_ind2, series2filt[shift:end] .- s2_avg)
   series2filt_ne = series2filt_slope*x_ind2 .+ s2_avg
   plot(series2time[shift:end], series2filt_ne, color = "r")
   xpos = series2time[1000]
   ypos = s2_avg
   xlabel("time")
   legend(["ERSST", "ERSST smoothed", "ERSST smoothed trend: " * string(round(series2filt_slope, digits =6))])
   suptitle(t)
   savefig(s)
   close()
end

"""
   function fftplot(ffte_freq, ffte, s, t)
   Plot frequency vs intensity for FFTW resuls
# Arguments#
- `ffte_freq`: output of FFTW.fft_freq
- `ffte`: ouptu of FFTW.fft
- `s`: save strng
- `t`:title string
"""
function fftplot(ffte_freq, ffte, s, t)
   load_python()
   figure()
   ax = subplot()
   N = length(ffte_freq)
   halfway = 0.0
   if N % 2 == 1
      halfway = convert(Int, N/2 - 0.5 )
   else
      halfway = convert(Int, N/2)
   end
   replace!(ffte_freq, 0.0 => 1.0)
   print(halfway)

   if typeof(ffte[1]) in [Vector{ComplexF32}, Vector{ComplexF64}]
      for i in ffte
         ax.plot(ffte_freq[begin+1:halfway], abs.(i)[begin+1:halfway], linewidth = 2, alpha = 0.5)
      end
      legend(["Raw", "Anomaly"])
   else
      ax.plot(ffte_freq[begin+1:halfway], abs.(ffte)[begin+1:halfway])
   end
   xt = range(ffte_freq[begin + 1], stop = ffte_freq[halfway], step = 2*10^(-10) * 10)
   ax.set_xticks(xt)
   ax.set_xlabel("Frequency [Hz]")
   ax.set_ylabel("Intensity")

   #ax.set_ylim([0, 750])
   ax.set_xlim(ffte_freq[begin+1], 2 * 0.1*10^(-7))
   ax.set_ylim(0, 500)
   ax.grid()
   title(t)
   savefig(s)
   close()

end

"""
   function normalequation(x,y)
   Linear regression with least squares fit
   Assumes y-intercept is already 0 (mean removed )
# Arguments#
- `x`: dependent variable
- `y`: independent variable
# Output
- `val`: slope of line of best fit
"""
function normalequation(x,y)
   val = inv(transpose(x)*x)*transpose(x)*y
   return val
end

"""
   function tempmean(array)
   Computes temporal mean of 3d array, IGNORING NANs
# Arguments#
- `array`: 3d float array in format (lat x lon x time)
# Output
- `avg`: 2d array in format (lat x lon), where every value is temporal average
"""
function temp_mean(array)
   avg = similar(array)[:,:,1]

   for i in range(1, stop = size(avg)[1], step = 1)
      for j in range(1, stop = size(avg)[2], step = 1)
         ts = array[i, j, :]
         avg[i, j] = NaNMath.mean(ts)
      end
   end
   return avg
end

"""
   function tempstd(array)
   Computes temporal standard deviation of 3d array, IGNORING NANs
# Arguments#
- `array`: 3d float array in format (lat x lon x time)
# Output
- `avg`: 2d array in format (lat x lon), where every value is temporal standard deviation
"""
function temp_std(array)
   avg = similar(array)[:,:,1]

   for i in range(1, stop = size(avg)[1], step = 1)
      for j in range(1, stop = size(avg)[2], step = 1)
         ts = array[i, j, :]
         avg[i, j] = NaNMath.std(ts)
      end
   end
   return avg
end

"""
   function spatial_mean(array)
   Computes spatial average of 3d array, IGNORING NANs
# Arguments#
- `array`: 3d float array in format (lat x lon x time)
# Output
- `avg`: 1d array in format (time), where every value is spatial average at that time
"""
function spatial_mean(array)
   avg = similar(array)[1,1,:]
   for i in range(1, stop = size(array)[3], step = 1)
      ts = array[:,:, i]
      ts = vcat(ts...)#flatten
      avg[i] = NaNMath.mean(ts)
   end
   return avg
end

"""
   function spatial_std(array)
   Computes spatial average of 3d array, IGNORING NANs
# Arguments#
- `array`: 3d float array in format (lat x lon x time)
# Output
- `avg`: 1d array in format (time), where every value is spatial standard deviation at that time
"""
function spatial_std(array)
   avg = similar(array)[1,1,:]
   for i in range(1, stop = size(array)[3], step = 1)
      ts = array[:,:, i]
      ts = vcat(ts...)#flatten
      avg[i] = NaNMath.std(ts)
   end
   return avg
end

function plot_eofs(solver, lon, lat, time, s)
   load_python()
   eofs = solver.eofs(neofs=2)
   pcs = solver.pcs(npcs = 2, pcscaling = 1)
   var = solver.varianceFraction()
   vmin = NaNMath.minimum(eofs)
   vmax = NaNMath.maximum(eofs)

   vmin = abs(vmin) > abs(vmax) ? vmin : -vmax
   vmax = abs(vmax) > abs(vmin) ? vmax : -vmin

   figure(figsize = (12, 5))
   ax = subplot(2,2,1, projection = ccrs.PlateCarree())
   gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=true,
                     linewidth=2, color="gray", alpha=0.5, linestyle="--")
   gl.right_labels = false
   gl.top_labels = false

   pc = contourf(lon, lat, transpose(eofs[1, :, :]), cmap = "coolwarm", transform = ccrs.PlateCarree(), vmin = vmin, vmax = vmax)
   frac = 0.013
   cb = colorbar(pc, fraction = frac)
   ax.coastlines()
   cb.set_label("SST Anomaly [σ]")
   title("EOF1, Exp. Var = " *string(round(var[1]*100, digits = 2)) * "%")
   ax = subplot(2,2,3)
   plot(time, pcs[:, 1])
   ax.set_ylim(-2.5, 2.5)
   xlabel("time [years]")
   ax.set_xlim(Dates.DateTime(1860,1,1), time[end])
   grid()
   title("PC1")

   ax = subplot(2,2,2, projection = ccrs.PlateCarree())
   gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=true,
                     linewidth=2, color="gray", alpha=0.5, linestyle="--")
   gl.right_labels = false
   gl.top_labels = false
   gl.left_labels = false
   pc = contourf(lon, lat, transpose(eofs[2, :, :]), cmap = "coolwarm", transform = ccrs.PlateCarree(), vmin = vmin, vmax = vmax )
   cb = colorbar(pc, fraction = frac)
   ax.coastlines()
   cb.set_label("SST Anomaly [σ]")
   title("EOF2, Exp. Var = " *string(round(var[2]*100, digits = 2)) * "%")
   ax = subplot(2,2,4)
   plot(time, pcs[:, 2])
   xlabel("time [years]")
   ax.set_ylim(-2.5, 2.5)
   ax.set_xlim(Dates.DateTime(1860,1,1), time[end])
   title("PC2")
   grid()
   savefig(s,bbox_inches="tight")
   close()
end

function plot_corr(xc, lat, lon, time, s)
   load_python()
   figure()
   ax = subplot(2,1,1, projection = ccrs.PlateCarree())
   ax.coastlines()
   halfway = (size(xc)[1] - 1)/2
   halfway = convert(Int32, halfway)
   vmin = NaNMath.minimum(xc)
   vmax = NaNMath.maximum(xc)

   vmin = abs(vmin) > abs(vmax) ? vmin : -vmax
   vmax = abs(vmax) > abs(vmin) ? vmax : -vmin

   pc = contourf(lon, lat, transpose(xc[halfway, :, :]), cmap = "coolwarm", transform = ccrs.PlateCarree(), vmin = vmin, vmax = vmax)
   cb = colorbar(pc, fraction = 0.013)
   cb.set_label("Correlation")
   gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=true,
                     linewidth=2, color="gray", alpha=0.5, linestyle="--")
   gl.right_labels = false
   gl.top_labels = false
   subplot(2,1,2)
   reshaped = zeros(size(xc)[2], size(xc)[3], size(xc)[1])

   permutedims!(reshaped, xc,[2,3,1])
   mean = spatial_mean(reshaped)
   lag = -halfway:halfway
   plot(lag, mean)
   xlabel("Lag [Months]")
   ylabel("Autocorrelation")
   savefig(s,bbox_inches="tight")
   close()
end

function plot_reg(reg, lat, lon, s)
   load_python()
   figure()
   ax = subplot(projection = ccrs.PlateCarree())
   ax.coastlines()
   vmin = NaNMath.minimum(reg)
   vmax = NaNMath.maximum(reg)

   vmin = vmin < 0 ? vmin : -vmin
   vmin = abs(vmin) > abs(vmax) ? vmin : -vmax
   vmax = abs(vmax) > abs(vmin) ? vmax : -vmin

   pc = contourf(lon, lat, transpose(reg), cmap = "coolwarm", vmin = vmin, vmax = vmax)
   cb = colorbar(pc, fraction = 0.013, format="%.0e")
   cb.set_label("Regression Coefficient")
   gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=true,
                     linewidth=2, color="gray", alpha=0.5, linestyle="--")
   gl.right_labels = false
   gl.top_labels = false

   savefig(s,bbox_inches="tight")
   close()
end

end


text:
"""
Flatten a vector of vectors
"""
function flatten{T}(x::Vector{Vector{T}})
    out = T[]
    for y in x
        append!(out, y)
    end
    out
end

"""
```
metisform(G) -> n, xadj, adjncy
```

Represent an undirected, unweighted graph `G` in form required for
Metis functions. `G` can be an adjacency matrix (sparse or dense), an
adjacency list from the Graphs package, or a graph from the the
LightGraphs package.  `n` is the number of vertices, `xadj` is an
array of 0-based column pointers, and `adjncy` is an array of 0-based
row indices.
"""
# Convert abstract matrix to Metis form
function metisform{T,N}(m::AbstractMatrix{N}, ::Type{T})
    metisform(convert(SparseMatrixCSC{N,T},m)n, T)
end

# Convert CSC matrix to METIS form
function metisform{T}(m::SparseMatrixCSC, ::Type{T})
    issymmetric(m) || ishermitian(m) || error("m must be symmetric or Hermitian")

    # Copy m.rowval and m.colptr to T vectors and drop diagonal elements
    adjncy = @compat sizehint!(T[],nnz(m))
    xadj = zeros(T,m.n+1)
    for j in 1:m.n
        count = 0
        for k in m.colptr[j] : (m.colptr[j+1] - 1)
            i = m.rowval[k]
            if i != j
                count += 1
                push!(adjncy,i-1)
            end
        end
        xadj[j+1] = xadj[j] + count
    end
    convert(T,m.n),xadj,adjncy
end

# Convert Graphs.GenericAdjacencyList to Metis form
function metisform{T}(al::Graphs.GenericAdjacencyList, ::Type{T})
    !Graphs.is_directed(al) || error("Metis functions require undirected graphs")
    @compat isa(al.vertices,UnitRange) && first(al.vertices) == 1 || error("Vertices must be numbered from 1")
    @compat length(al.adjlist), round(T, cumsum(vcat(0, map(length, al.adjlist)))),
        round(T, flatten(al.adjlist)) .- one(T)
end

# Convert LightGraphs.Graph to Metis form
function metisform{T}(g::LightGraphs.Graph, ::Type{T})
    n = convert(T,nv(g))
    xadj = zeros(T,n + 1)
    adjncy = sizehint!(T[],2*ne(g))
    for i in 1:n
        ein = [convert(T,dst(x)-1) for x in LightGraphs.out_edges(g, i)]
        xadj[i + 1] = xadj[i] + length(ein)
        append!(adjncy,ein)
    end
    n, xadj, adjncy
end

"""
```
metisform_weighted(G) -> n, xadj, adjncy, adjwgt
```

Represent an undirected, weighted graph `G` in form required for Metis
functions. `G` should be an adjacency matrix (sparse or dense).  `n`
is the number of vertices, `xadj` is an array of 0-based column
pointers, `adjncy` is an array of 0-based row indices, and `adjwgt` is
an array of edge weights. Edge weights are rounded to the nearest
integer since Metis requires integer edge weights.
"""
# Convert abstract matrix to Metis form
function metisform_weighted{T,N}(m::AbstractMatrix{N}, ::Type{T})
    metisform_weighted(convert(SparseMatrixCSC{N,T},m), T)
end

# Convert CSC matrix to Metis form
function metisform_weighted{T}(m::SparseMatrixCSC, ::Type{T})
    issymmetric(m) || ishermitian(m) || error("m must be symmetric or Hermitian")

    # copy m.rowvalm, m.colptr, and m.nzval to T vectors
    # and drop diagonal elements
    adjncy = @compat sizehint!(T[],nnz(m))
    adjwgt = @compat sizehint!(T[],nnz(m))
    xadj = zeros(T,m.n+1)
    for j in 1:m.n
        count = 0
        for k in m.colptr[j] : (m.colptr[j+1] - 1)
            i = m.rowval[k]
            if i != j
                count += 1
                push!(adjncy,i-1)
                push!(adjwgt,round(T,m[i,j]))
            end
        end
        xadj[j+1] = xadj[j] + count
    end
    convert(T,m.n),xadj,adjncy,adjwgt
end

"""
```
testgraph(nm) -> Graphs.GenericAdjacencyList
```

Load a file in the `Metis.jl/graphs/` directory and construct a
graph. `nm` should be the file name, excluding the `.graph` file
extension.
"""
# Load test graph
function testgraph(nm::String)
    pathnm = joinpath(dirname(@__FILE__), "..", "graphs", string(nm, ".graph"))
    ff = open(pathnm, "r")
    nvert, nedge = map(t -> parse(Int, t), split(readline(ff)))
    adjlist = Array(Vector{Int32}, nvert)
    for i in 1:nvert adjlist[i] = map(t -> parse(Int32, t), split(readline(ff))) end
    close(ff)
    @compat GenericAdjacencyList{Int32,UnitRange{Int32},Vector{Vector{Int32}}}(false,
                                                                    map(Int32, 1:nvert),
                                                                    nedge,
                                                                    adjlist)
end


text:
module IDXParser

# package code goes here

using Base.Docs

export parseIDX

@doc """
Parse IDX-formatted file and return an array of values:
  1ˢᵗ byte: 0x00
  2ⁿᵈ byte: 0x00
  3ʳᵈ byte: Array Datatype: 0x08 -> UInt8, 0x09 -> Int8, 0x0B -> Int16, 0x0C -> Int32, 0x0D -> Float32, 0x0E -> Float64
  4ᵗʰ byte: # dimensions: 0x01 for vector, 0x02 for 2D matrix, etc.
  5ᵗʰ byte: length of first dimension
  ⋮
  (5 + #dimensions - 1)ᵗʰ byte: length of final dimension
  (5 + #dimensions)ᵗʰ byte: start of data
""" ->
function parseIDX(file_name)
    fs = open(file_name, "r")
    magic_number = read(fs, UInt8, 4)

    if magic_number[1] != 0x00 && magic_number[2] != 0x00
        throw(DomainError())
    end

    if magic_number[3] == 0x08
        idx_type = UInt8
    elseif magic_number[3] == 0x09
        idx_type = Int8
    elseif magic_number[3] == 0x0B
        idx_type = Int16
    elseif magic_number[3] == 0x0C
        idx_type = Int32
    elseif magic_number[3] == 0x0D
        idx_type = Float32
    elseif magic_number[3] == 0x0E
        idx_type = Float64
    else
        throw(DomainError())
    end

    dims = Array{Int64}(map(ntoh, read(fs, Int32, magic_number[4])))

    contents = read(fs, idx_type, prod(dims))

    if idx_type != UInt8 && idx_type != Int8
        contents = map(ntoh, contents)
    end

    out = reshape(contents, tuple(dims...))
    return out
end

end # module


text:
"""
	@author: Alois Pichler
	Created September 2018
		provides ontoSimplex: function to project onto the simplex
"""

function ontoSimplex!(probabilities::Vector{Float64})   # projects the probabilities onto the simplex
    for i= 1:length(probabilities)
        if probabilities[i] < 0.0   # make them nonnegative
            probabilities[i]= 0.0
        end
    end
	summ= sum(probabilities)
	if summ <= 0.0
		probabilities= fill!(probabilities, 1.0/ length(probabilities))
		summ= sum(probabilities)
		@info "ontoSimplex: created probabilities." maxlog=1
	end
	if summ != 1
		probabilities./= summ
		tmpi= rand(1:length(probabilities)) 	# modify a random index
		summ= -sum(probabilities[1:tmpi-1])+ 1.0- sum(probabilities[tmpi+1:end])
		probabilities[tmpi]= max(0, summ)
		@info "ontoSimplex: modified probabilities." maxlog=1
	end
    nothing     # modifies probabilities
end


text:
##LossFunction type: define the loss function used and its hyper-parameter
@compat abstract type LossFunction end
@compat abstract type Regression <: LossFunction end
@compat abstract type Classification <: LossFunction end
  #Loss functions for regression
  struct OLS <: Regression
  end
  struct L1SVR <: Regression
    ɛ::Float64
  end
  struct L2SVR <: Regression
    ɛ::Float64
  end

  #Loss functions for classification
  struct LogReg <: Classification
  end
  struct L1SVM <: Classification
  end
  struct L2SVM <: Classification
  end

  ##Point-wise value of the Fenchel conjugate for each loss function
  function loss(ℓ::OLS, y, u)
    return .5*(y.-u).^2
  end
  function loss(ℓ::L1SVR, y, u)
    return max.(abs.(y-u) .- ℓ.ɛ, 0.)
  end
  function loss(ℓ::L2SVR, y, u)
    return .5*max.(abs.(y-u) .- ℓ.ɛ, 0.)^2
  end
  function loss(ℓ::LogReg, y, u)
    return log.(1 .+ exp.(-y.*u))
  end
  function loss(ℓ::L1SVM, y, u)
    return max.(0., 1 .- y.*u)
  end
  function loss(ℓ::L2SVM, y, u)
    return .5*max.(0., 1 .- y.*u).^2
  end

  ##Point-wise value of the Fenchel conjugate for each loss function
  function fenchel(ℓ::OLS, y, a)
    return .5*a^2 + a*y
  end
  function fenchel(ℓ::L1SVR, y, a)
    return y*a + ℓ.ɛ*abs(a)
  end
  function fenchel(ℓ::L2SVR, y, a)
    return .5*a^2 + y*a + ℓ.ɛ*abs(a)
  end
  function fenchel(ℓ::LogReg, y, a)
    return (1+a*y)*log(1+a*y) - a*y*log(-a*y)
  end
  function fenchel(ℓ::L1SVM, y, a)
    return a*y
  end
  function fenchel(ℓ::L2SVM, y, a)
    return .5*a^2 + a*y
  end

  ##Point-wise derivative of the Fenchel conjugate for each loss function
  function grad_fenchel(ℓ::OLS, y, a)
    return a + y
  end
  function grad_fenchel(ℓ::L1SVR, y, a)
    return y + ℓ.ɛ*sign(a)
  end
  function grad_fenchel(ℓ::L2SVR, y, a)
    return a + y + ℓ.ɛ*sign(a)
  end
  function grad_fenchel(ℓ::LogReg, y, a)
    return y*log(1+a*y) - y*log(-a*y)
  end
  function grad_fenchel(ℓ::L1SVM, y, a)
    return y
  end
  function grad_fenchel(ℓ::L2SVM, y, a)
    return a + y
  end

##Sparsity type: specify how sparsity is enforced, constrained or penalized
@compat abstract type Sparsity end
    #Constraint: add the constraint "s.t. ||w||_0<=k"
    struct Constraint <: Sparsity
      k::Int
    end
    function parameter(Card::Constraint)
      return Card.k
    end
    max_index_size(Card::Constraint, p::Int) = min.(Card.k, p)

    #Penalty: add the penalty "+λ ||w||_0"
    struct Penalty <: Sparsity
      λ::Float64
    end
    function parameter(Card::Penalty)
      return Card.λ
    end
    max_index_size(Card::Penalty, p::Int) = p

##SparseEstimator type: output of the algorithm
struct SparseEstimator
  "Loss function used in the model"
  loss::LossFunction
  "Model to account for sparsity"
  sparsity::Sparsity
  "Set of relevant indices"
  indices::Array
  "Estimator w on those selected indices"
  w::Array{Float64}
  "Dual variables α"
  α::Array{Float64}
  "Intercept/bias term"
  b::Float64
  "Number of iterations in the sug-gradient algorithm"
  iter::Integer
end


text:
using OpenTelemetryAPI
using OpenTelemetrySDK
using OpenTelemetryProto
using OpenTelemetryExporterOtlpProtoGrpc

using Documenter

open(joinpath(@__DIR__, "..", "README.md"), "r") do f_src
    open(joinpath(@__DIR__, "src", "index.md"), "w") do f_dest
        s_dest = read(f_src, String)
        s_dest = replace(s_dest, "<!-- ```@raw html -->" => "```@raw html")
        s_dest = replace(s_dest, "<!-- ``` -->" => "```")
        s_dest = replace(s_dest, "```julia" => "```@example")
        write(f_dest, s_dest)
    end
end

makedocs(
    modules = [
        OpenTelemetryAPI
        OpenTelemetrySDK
        OpenTelemetryProto
        OpenTelemetryExporterOtlpProtoGrpc
    ],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "G-YWC6SQHSVM",
        assets = ["assets/favicon.ico"],
    ),
    sitename = "OpenTelemetry.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Design" => Any["API"=>"design_api.md", "SDK"=>"design_sdk.md"],
    ],
)

deploydocs(repo = "github.com/oolong-dev/OpenTelemetry.jl.git", target = "build")


text:
## PART 0

# Re initialize stuff from question 1

println("\n~~ QUESTION 2 ~~\n");
println("RUNNING QUESTION 1 SCRIPT AGAIN");

# Load in the datasets and store in array v
include("../q1/p1.jl");

# Take means of the datasets and store in array v̄1 and v̄2
include("../q1/p2.jl");


text:
# various recipes for Plots.jl

"""
   trimplot(f, a, b, c=20; kwargs...)

Plot f over [a,b] but break graph if it exceeds c in absolute value.
"""
@userplot TrimPlot
@recipe function __(a::TrimPlot)
    if length(a.args) < 4
        f, x0, x1 = a.args
        c = 20
    else
        f, x0, x1, c = a.args
    end

    xs = range(x0, x1, length=251)
    ys = f.(xs)

    Xs = Any[]
    Ys = Any[]

    us, vs = Real[], Real[]
    for (x,y) in zip(xs, ys)
        if abs(y) <= c
            push!(us, x); push!(vs, y)
        else
            if length(us) > 0
                push!(Xs, copy(us))
                push!(Ys, copy(vs))
            end
            empty!(us); empty!(vs)
        end
    end
    if length(us) > 0
        push!(Xs, copy(us))
        push!(Ys, copy(vs))
    end

    seriestype := :path
    x := Xs
    y := Ys

    Xs, Ys
end

"""
    plotif(f, g, a, b)

Plot f colored depending on g < 0 or not.
"""
@userplot PlotIf
@recipe function __(a::PlotIf)
    f, g, x0, x1 = a.args
    xs = range(x0, x1, length=251)
    ys = f.(xs)

    seriestype := :path
    x := xs
    y := ys
    color := [gx < 0 ? :red : :blue for gx in g.(xs)]
    linewidth := 5
    legend := false

    xs, ys
end

@userplot SignChart
@recipe function __(a::SignChart)
    f, x0, x1 = a.args
    xs = range(x0, x1, length=251)
    ys = f.(xs)
    zs = zeros(length(xs))
    @series begin
        seriestype := :path
        x := xs
        y := zs
        linewidth := 1
        color := :black
        xs, zs
    end

    seriestype := :path
    x := xs
    y := ys
    color := [y < 0 ? :red : :blue for y in ys]

    linewidth --> 5
    legend --> false

    xs, ys
end

#### Multidimensional

# """
#     plot_polar(r, a, b; [n=251]

# Create a polar plot of `r(theta)` for `a <= theta <= b`.

# An alternative to the three step approach:

# ```
# ts = range(a, b, length=n)
# rs = r.(ts)
# plot(ts, rs, proj=:polar, legend=false)
# ```

# """
# @userplot Plot_Polar
# @recipe function __(rt::Plot_Polar; n=251)

#     r, a, b = rt.args
#     fn = t -> r(t) * [cos(t), sin(t)]
#     ts = range(a, stop=b, length=n)
#     xs = fn.(ts)

#     xys = unzip(xs)
#     seriestype := :path
#     x := xys[1]
#     y := xys[2]
#     xys
# end

# """
#     plot_parametric_curve(r, a, b; [n=251])

# For `r` taking `R -> R^n` make a plot over [a,b] using n points.

# There are other ways to produce this plot:
# * if we have two functions, `rx` and `ry` say, then passing as function arguments should work

# ```
# plot(rx, ry, a, b)
# ```

# * if we have a vector-valued function, then either of these work:

# ```
# r(t) = [cos(t), sin(t), t]
# plot(unzip(r, a, b, n)...)  # unzip reshapes the data
# plot_parametric_curve(r, a, b)
# ```

# We see that this function is a convenience for vector-valued descriptions of parameterized curves.
# """
# @userplot Plot_Parametric_Curve
# @recipe function __(r::Plot_Parametric_Curve; n=251)

#     f, a, b = r.args
#     ts = range(a, stop=b, length=n)
#     xs = f.(ts)

#     n = length(first(xs))
#     if n == 2
#         xys = unzip(xs)
#         seriestype := :path
#         x := xys[1]
#         y := xys[2]
#         xys
#     else
#         xyzs = unzip(xs)
#         seriestype := :path3d
#         x := xyzs[1]
#         y := xyzs[2]
#         z := xyzs[3]
#         xyzs
#     end
# end

# """
#     plot_parametric_surface(F, [xlim=(-5,5)], [ylim=(-5,5)], [nx=50], [ny=50])

# Plot the paremeterized surface described by `F:R^2 -> R^3`. This works with `pyplot` and `plotly`, but not `gr`.

# A surface `(x,y,f(x,y))` can be directly plotted with `surface` as:

# ```
# surface(xs, ys, f)
# ```

# where `xs` and `ys` provide a grid to plot over.

# For parametrically described surfaces, the above doesn't work. This function provides an interface:

# ```
# F(u,v) = [u*cos(v), u*sin(v), u] # a cone
# plot_parametric_surface(F, xlim=(0, 1), ylim=(0,2pi))
# ```

# The values of `xlim` are used as the range to plot over; `nx` specifies the number of points to use. Similarly for `ylim`.
# """
# @userplot Plot_Parametric_Surface
# @recipe function __(r::Plot_Parametric_Surface; nx=50, ny=50)

#     F = first(r.args)

#     xlim = get(plotattributes,:xlims, (-5,5))
#     ylim = get(plotattributes,:ylims, (-5,5))

#     us = range(xlim[1], xlim[2], length=nx)
#     vs = range(ylim[1], ylim[2], length=ny)
#     ws = unzip(F.(us, vs'))

#     seriestype := :surface
#     x := ws[1]
#     y := ws[2]
#     z := ws[3]
#     xlims := extrema(ws[1])
#     ylims := extrema(ws[2])
#     ws

# end

"""
   arrow(p, v)

Add vector, `v`, to plot anchored at point `p`.

Example
```
Fn = parametric(t -> [2cos(t), 3sin(t)])
Fnp = t -> ForwardDiff.derivative(Fn, t)
p = plot(Fn, 0, 2pi, legend=false)
for t in 0:pi/4:pi
   arrow!(Fn(t), Fnp(t))
end
p
```
"""
@userplot Arrow
@recipe function __(a::Arrow)
    v0, dv0 = a.args
    v = [v0...]; dv = [dv0...]
    if length(v) == 2
        seriestype := :quiver
        quiver := unzip([dv])
        x,y = unzip([v])
        x := x
        y := y
        ()
    else
        seriestype := :path3d
        M = hcat(v, v + dv)'
        (M[:,1],M[:,2], M[:,3])
    end

end

"""
    vectorfieldplot(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=8], [ny=8])

Create a vector field plot using a grid described by `xlim`, `ylim` with `nx` and `ny` grid points in each direction.

```
F(x,y) = [-y, x]
vectorfieldplot(F, xlim=(-4,4), ylim=(-4,4))
```
"""
@userplot VectorFieldPlot
@recipe function __(V::VectorFieldPlot; nx=8, ny=8)

    ars = V.args
    F = ars[1]
    xlim = get(plotattributes,:xlims, (-5,5))
    ylim = get(plotattributes,:ylims, (-5,5))

    dx, dy = (xlim[2]-xlim[1])/(nx-1), (ylim[2]-ylim[1])/(ny-1)
    xs, ys = xlim[1]:dx:xlim[2], ylim[1]:dy:ylim[2]

    ps = [[x,y] for x in xs for y in ys]
    vs = (xy -> F(xy[1], xy[2])).(ps)

    λ = 0.9 * minimum([u/maximum(getindex.(vs,i)) for (i,u) in enumerate((dx,dy))])

    seriestype := :quiver
    quiver := unzip(λ * vs)
    xs, ys = unzip(ps)
        x := xs
    y := ys
    ()
end

"""
    vectorfieldplot3d(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=5], [ny=5])

Create a 3 dimensional vector field plot using a grid described by `xlim`, `ylim`, `zlim` with `nx`, `ny`, and `nz` grid points in each direction.

Note: the vectors are represented with line, not arrow due to no implementation of `:quiver3d`.

```
F(x,y,z) = [-y, x,z]
vectorfieldplot3d(F, xlims=(-4,4), ylims=(-4,4), zlims=(0,3))
```
"""
@userplot VectorFieldPlot3d
@recipe function __(V::VectorFieldPlot3d;  nx=5, ny=5, nz=5)
    ars = V.args
    F = ars[1]
    xlim = get(plotattributes,:xlims, (-5,5))
    ylim = get(plotattributes,:ylims, (-5,5))
    zlim = get(plotattributes,:zlims, (-5,5))

    dx, dy, dz = (xlim[2]-xlim[1])/(nx-1), (ylim[2]-ylim[1])/(ny-1), (zlim[2]-zlim[1])/(nz-1)
    xs, ys, zs = xlim[1]:dx:xlim[2], ylim[1]:dy:ylim[2], zlim[1]:dz:zlim[2]

    ps = [[x,y,z] for x in xs for y in ys for z in zs]
    vs = (xyz -> F(xyz[1], xyz[2], xyz[3])).(ps)

    λ = 0.9 * minimum([u/maximum(getindex.(vs,i)) for (i,u) in enumerate((dx,dy,dx))])

    seriestype := :path3d
    dxs, dys, dzs = unzip(λ * vs)
    xs, ys, zs = unzip(ps)

    legend := false
    [[x,dx] for (x,dx) in zip(xs, xs+dxs)],    [[y,dy] for (y,dy) in zip(ys, ys+dys)],    [[z,dz] for (z,dz) in zip(zs, zs+dzs)]
end


text:
module PGField

export GFElem, PGInf, order

struct GFElem{p}
    value::Int
    GFElem{p}(value::Integer) where p = new{p}(mod(value, p))
end
struct PGInf{p} end
Base.zero(::Union{GFElem{p},PGInf{p}}) where p = GFElem{p}(0)
Base.one(::Union{GFElem{p},PGInf{p}}) where p = GFElem{p}(1)
Base.iszero(x::GFElem) = iszero(x.value)
Base.iszero(::PGInf) = false

# Addition
Base.:+(x::GFElem{p},y::GFElem{p}) where p = GFElem{p}(widen(x.value)+widen(y.value))
Base.:+(x::GFElem{p},∞::PGInf{p}) where p = ∞
Base.:+(∞::PGInf{p},x::GFElem{p}) where p = x+∞
Base.:+(∞::PGInf{p},::PGInf{p}) where p = zero(∞)

# Multiplication
Base.:*(x::GFElem{p},y::GFElem{p}) where p = GFElem{p}(widemul(x.value,y.value))
Base.:*(x::GFElem{p},∞::PGInf{p}) where p = iszero(x) ? one(x) : ∞
Base.:*(∞::PGInf{p},x::GFElem{p}) where p = x*∞
Base.:*(∞::PGInf{p},::PGInf{p}) where p = ∞  # is this right?

# Negation and inverse
Base.:-(x::GFElem{p}) where p = GFElem{p}(-x.value)
Base.:-(∞::PGInf) = ∞
function Base.inv(x::GFElem{p}) where p  # based on https://github.com/Nemocas/AbstractAlgebra.jl/blob/master/src/julia/GF.jl#L270-L277
    iszero(x) && return PGInf{p}()
    g, s, t = gcdx(x.value, p)
    g != 1 && error("Characteristic not prime in GF($p).")
    return GFElem{p}(s)
end
Base.inv(∞::PGInf) = zero(∞)

# Powers
Base.:^(x::GFElem{p},n) where p = GFElem{p}(BigInt(x.value)^n)

# Order of element
order(x::GFElem{p}) where p = length(unique(x^i for i in 1:p))

# Broadcasting
Base.broadcastable(x::GFElem) = Ref(x)
Base.broadcastable(x::PGInf) = Ref(x)

end # module


text:
# The PotentialFlow.jl package is licensed under the MIT "Expat" License:
#
# > Copyright (c) 2017: Darwin Darakananda.
# >
# > Permission is hereby granted, free of charge, to any person obtaining a copy
# > of this software and associated documentation files (the "Software"), to deal
# > in the Software without restriction, including without limitation the rights
# > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# > copies of the Software, and to permit persons to whom the Software is
# > furnished to do so, subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be included in all
# > copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.
# >

# """
# A macro for extracting fields from an object.  For example, instead of a statement
# like
#
#     (obj.a + obj.b)^obj.c
#
# it is more readable if we first locally bind the variables
#
#     a = obj.a
#     b = obj.b
#     c = obj.c
#
#     (a + b)^c
#
# Using the `@get` marco, this becomes
#
#     @get obj (a, b, c)
#
#     (a + b)^c
#
# We can also locally assign different names to the binding
#
#     @get obj (a, b, c) (α, β, γ)
#     (α + β)^γ
# """
macro get(object, fields...)
    if length(fields) == 1
        try
            @assert typeof(fields[1]) == Expr
            @assert fields[1].head == :tuple
            @assert all([typeof(arg) == Symbol for arg in fields[1].args])
        catch
            throw(ArgumentError("second argument must be a tuple of field names"))
        end
        esc(Expr(:block, [:($sym = $object.$sym) for sym in fields[1].args]...))
    elseif length(fields) == 2
        try
            @assert typeof(fields[1]) == Expr
            @assert typeof(fields[2]) == Expr
            @assert (fields[1].head == :tuple && fields[2].head == :tuple)
            @assert all([typeof(arg) == Symbol for arg in fields[1].args])
            @assert all([typeof(arg) == Symbol for arg in fields[2].args])
        catch
            throw(ArgumentError("second and third argument must be tuples of field names"))
        end

        nargs = length(fields[1].args)
        @assert nargs == length(fields[2].args) "field name tuples must have the same length"
        esc(Expr(:block, [:($(fields[2].args[i]) = $object.$(fields[1].args[i])) for i in 1:nargs]...))
    else
        throw(ArgumentError("Usage: @get <object> (field names...) [(reference names...)]"))
    end
end


text:
using ApproxFun

###
# Fractional Derivative
###

S = Legendre()

n = 50; # number of data points
p = range(-1,stop=1,length=n);   # a non-default grid

v = exp.(p);           # values at the non-default grid
m = 20   # number of basis functions
V = Array{Float64}(undef,n,m); # Create a Vandermonde matrix by evaluating the basis at the grid
for k = 1:size(V,2)
   V[:,k] = Fun(S,[zeros(k-1);1]).(p)
end
f = Fun(S,V\v); # least squares fit to data

f_rl = LeftDerivative(0.5)*f # RL left half-derivative
f_cap = LeftIntegral(0.5)*Fun(f', S) # Caputo left half-derivative

using Plots
plot(f_rl;ylims=(0,3),legend=:bottomright,label="RL")
plot!(f_cap; label="Caputo")
savefig("RLvCaputo.png")

max_stars_repo_path: Extras/Fractional derivatives.jl

max_issues_repo_path: Extras/Fractional derivatives.jl

max_forks_repo_path: Extras/Fractional derivatives.jl