text:
"""
This file stores the analysis functions of ClimateSatellite.jl that are important for
climatological analysis of the downloaded data.

Current functionalities stored in this file include:
- Yearly, seasonal and monthly means, maximum, minimum, range and standard deviation:
    - On a general climatological basis
    - For a specified hour
    - Of the diurnal variability
- Zonal and meridional averages of the above
- Saving the results of the abovementioned analysis into yearly NetCDF files
- Retrieval of these data from their respective NetCDF files

Functionalities that are in development include:
- Zonal, meridional and domain averages

Because there is progressively more data every year, analysis will be conducted on a yearly
basis instead of for the entire time-series.  This is to ensure generality and ease-of-use
as the datasets become larger and larger.

"""

function clisatopNaN(f::Function,data::Vector{Int16})
    ndata = @view data[data.!=-32768];
    if ndata != []; return f(ndata); else; return -32768; end
end

function clisatanalysis(
    productID::AbstractString, yr::Integer;
    varname::AbstractString,
    path::AbstractString="", region::AbstractString="GLB"
)

    if path == ""; dataroot = clisatroot(productID);
    else;          dataroot = clisatroot(productID,path);
    end
    isgeoregion(region);

    info = Dict{Any,Any}("root"=>dataroot); clisatinfo!(info,productID);
    clisatvarinfo!(info,productID,varname=varname);

    regfol = clisatregfol(info,region);
    rawfol = clisatrawfol(info,Date(yr),region);
    cd(rawfol);

    @info "$(Dates.now()) - Extracting $(info["source"]) $(info["product"]) data in $(gregionfullname(region)) region during $yr ..."

    lon,lat = clisatlonlat(info); rlon,rlat,rinfo = gregiongridvec(region,lon,lat);
    nlon = length(rlon); nlat = length(rlat); nt = info["dayfreq"]+1; grid = [rlon,rlat];

    davg = rand(Int16,nlon,nlat,nt+1,13); dstd = rand(Int16,nlon,nlat,nt+1,13);
    dmax = rand(Int16,nlon,nlat,nt+1,13); dmin = rand(Int16,nlon,nlat,nt+1,13);

    zavg = rand(Int16,nlat,nt+1,13); zstd = rand(Int16,nlat,nt+1,13);
    zmax = rand(Int16,nlat,nt+1,13); zmin = rand(Int16,nlat,nt+1,13);

    mavg = rand(Int16,nlon,nt+1,13); mstd = rand(Int16,nlon,nt+1,13);
    mmax = rand(Int16,nlon,nt+1,13); mmin = rand(Int16,nlon,nt+1,13);

    for mo = 1 : 12; ndy = daysinmonth(yr,mo)

        @info "$(Dates.now()) - Analyzing $(info["source"]) $(info["product"]) data in $(gregionfullname(region)) during $(Dates.monthname(mo)) $yr ..."
        ncraw = clisatrawname(productID,Date(yr,mo),region);
        ds = Dataset(ncraw,"r"); vds = ds[varname];
        raw = reshape(vds.var[:],nlon,nlat,(nt-1),ndy);
        close(ds);

        @debug "$(Dates.now()) - Extracting hourly information for each month ..."
        for it = 1 : nt-1, ilat = 1 : nlat, ilon = 1 : nlon

            rawii = @view raw[ilon,ilat,it,:]; rawii = @view rawii[rawii.!=-32768];

            if rawii != [];
                davg[ilon,ilat,it,mo] = round(Int16,mean(rawii));
                dstd[ilon,ilat,it,mo] = round(Int16,std(rawii));
                dmax[ilon,ilat,it,mo] = round(Int16,maximum(rawii));
                dmin[ilon,ilat,it,mo] = round(Int16,minimum(rawii));
            else
                davg[ilon,ilat,it,mo] = -32768; dstd[ilon,ilat,it,mo] = -32768;
                dmax[ilon,ilat,it,mo] = -32768; dmin[ilon,ilat,it,mo] = -32768;
            end

        end

        @debug "$(Dates.now()) - Permuting days and hours dimensions ..."
        raw = permutedims(raw,(1,2,4,3)); drg = zeros(nlon,nlat,ndy);
        dmn = zeros(nlon,nlat,ndy);

        @debug "$(Dates.now()) - Averaging diurnal data into daily data ..."
        for idy = 1 : ndy, ilat = 1 : nlat, ilon = 1 : nlon

            tmpii = @view raw[ilon,ilat,idy,:]; tmpii = @view tmpii[tmpii.!=-32768];

            if tmpii != [];
                drg[ilon,ilat,idy] = maximum(tmpii)/2 - minimum(tmpii)/2;
                dmn[ilon,ilat,idy] = mean(tmpii);
            else
                drg[ilon,ilat,idy] = -32768;
                dmn[ilon,ilat,idy] = -32768;
            end

        end

        @debug "$(Dates.now()) - Extracting information on monthly climatology and diurnal variability ..."
        for ilat = 1 : nlat, ilon = 1 : nlon;

            dmnii = @view dmn[ilon,ilat,:]; dmnii = @view dmnii[dmnii.!=-32768];
            drgii = @view drg[ilon,ilat,:]; drgii = @view drgii[drgii.!=-32768];

            davg[ilon,ilat,nt,mo] = round.(Int16,mean(dmnii));
            dstd[ilon,ilat,nt,mo] = round.(Int16,std(dmnii));
            dmax[ilon,ilat,nt,mo] = round.(Int16,maximum(dmnii));
            dmin[ilon,ilat,nt,mo] = round.(Int16,minimum(dmnii));

            davg[ilon,ilat,nt+1,mo] = round.(Int16,mean(drgii));
            dstd[ilon,ilat,nt+1,mo] = round.(Int16,std(drgii));
            dmax[ilon,ilat,nt+1,mo] = round.(Int16,maximum(drgii));
            dmin[ilon,ilat,nt+1,mo] = round.(Int16,minimum(drgii));

        end

    end

    @info "$(Dates.now()) - Calculating yearly climatology for $(info["source"]) $(info["product"]) in $(gregionfullname(region)) during $yr ..."
    for it = 1 : nt+1, ilat = 1 : nlat, ilon = 1 : nlon

        davg[ilon,ilat,it,end] = round.(Int16,clisatopNaN(mean,davg[ilon,ilat,it,1:12]));
        dstd[ilon,ilat,it,end] = round.(Int16,clisatopNaN(mean,dstd[ilon,ilat,it,1:12]));
        dmax[ilon,ilat,it,end] = round.(Int16,clisatopNaN(mean,dmax[ilon,ilat,it,1:12]));
        dmin[ilon,ilat,it,end] = round.(Int16,clisatopNaN(mean,dmin[ilon,ilat,it,1:12]));

    end

    @info "$(Dates.now()) - Calculating zonal-averaged climatology for $(info["source"]) $(info["product"]) in $(gregionfullname(region)) during $yr ..."
    for ilat = 1 : nlat, it = 1 : nt+1, imo = 1 : 13
        zavg[ilat,it,imo] = round.(Int16,clisatopNaN(mean,davg[:,ilat,it,imo]));
        zstd[ilat,it,imo] = round.(Int16,clisatopNaN(mean,dstd[:,ilat,it,imo]));
        zmax[ilat,it,imo] = round.(Int16,clisatopNaN(mean,dmax[:,ilat,it,imo]));
        zmin[ilat,it,imo] = round.(Int16,clisatopNaN(mean,dmin[:,ilat,it,imo]));
    end

    @info "$(Dates.now()) - Calculating meridional-averaged climatology for $(info["source"]) $(info["product"]) in $(gregionfullname(region)) during $yr ..."
    for imo = 1 : 13, it = 1 : nt+1, ilon = 1 : nlon
        mavg[ilon,it,imo] = round.(Int16,clisatopNaN(mean,davg[ilon,:,it,imo]));
        mstd[ilon,it,imo] = round.(Int16,clisatopNaN(mean,dstd[ilon,:,it,imo]));
        mmax[ilon,it,imo] = round.(Int16,clisatopNaN(mean,dmax[ilon,:,it,imo]));
        mmin[ilon,it,imo] = round.(Int16,clisatopNaN(mean,dmin[ilon,:,it,imo]));
    end

    clisatanasave([davg,dstd,dmax,dmin],[zavg,zstd,zmax,zmin],[mavg,mstd,mmax,mmin],
                  grid,productID,varname,yr,path,region,info);

end

function clisatanasave(
    data::Array{Array{Int16,4},1},
    zdata::Array{Array{Int16,3},1},
    mdata::Array{Array{Int16,3},1},
    grid::Vector{<:Any},
    productID::AbstractString, varname::AbstractString, yr::Integer,
    path::AbstractString, region::AbstractString, info::Dict
)

    @info "$(Dates.now()) - Saving analysed $(info["source"]) $(info["product"]) data in $(gregionfullname(region)) for the year $yr ..."

    fol = clisatanafol(info,region)
    fnc = joinpath(fol,clisatananame(productID,varname,Date(yr),region));
    rlon,rlat = grid; nlon = length(rlon); nlat = length(rlat); nt = info["dayfreq"];

    if isfile(fnc)
        @info "$(Dates.now()) - Stale NetCDF file $(fnc) detected.  Overwriting ..."
        rm(fnc);
    end

    @debug "$(Dates.now()) - Creating NetCDF file $(fnc) for analyzed $(info["source"]) $(info["product"]) data in $yr ..."

    ds = Dataset(fnc,"c");
    ds.dim["longitude"] = nlon; ds.dim["latitude"] = nlat;
    ds.dim["hour"] = nt; ds.dim["month"] = 12;

    att_var = map(x->Dict(),1:2)
    for ii = 1 : 2
        att_var[ii]["units"]         = info["units"];
        att_var[ii]["standard_name"] = info["standard"];
        att_var[ii]["long_name"]     = info["variable"];
        att_var[ii]["missing_value"] = -32768;
    end
    att_var[1]["scale_factor"] = info["scale"];   att_var[1]["add_offset"] = info["offset"];
    att_var[2]["scale_factor"] = info["scale"]*2; att_var[2]["add_offset"] = 0;

    att_lon = Dict("units"=>"degrees_east","long_name"=>"longitude");
    att_lat = Dict("units"=>"degrees_north","long_name"=>"latitude");

    dlon = defVar(ds,"longitude",Float32,("longitude",),attrib=att_lon); dlon[:] = rlon;
    dlat = defVar(ds,"latitude",Float32,("latitude",),attrib=att_lat);   dlat[:] = rlat;

    @debug "$(Dates.now()) - Saving analyzed $(info["source"]) $(info["product"]) data for $yr to NetCDF file $(fnc) ..."

    v = defVar(ds,"domain_yearly_mean_climatology",Int16,
               ("longitude","latitude"),attrib=att_var[1]);
    v.var[:] = data[1][:,:,nt+1,end];

    v = defVar(ds,"domain_yearly_std_climatology",Int16,
               ("longitude","latitude"),attrib=att_var[1]);
    v.var[:] = data[2][:,:,nt+1,end];

    v = defVar(ds,"domain_yearly_maximum_climatology",Int16,
               ("longitude","latitude"),attrib=att_var[1]);
    v.var[:] = data[3][:,:,nt+1,end];

    v = defVar(ds,"domain_yearly_minimum_climatology",Int16,
               ("longitude","latitude"),attrib=att_var[1]);
    v.var[:] = data[4][:,:,nt+1,end];

    v = defVar(ds,"domain_yearly_mean_hourly",Int16,
               ("longitude","latitude","hour"),attrib=att_var[1]);
    v.var[:] = data[1][:,:,1:nt,end];

    v = defVar(ds,"domain_yearly_std_hourly",Int16,
           ("longitude","latitude","hour"),attrib=att_var[1]);
    v.var[:] = data[2][:,:,1:nt,end];

    v = defVar(ds,"domain_yearly_maximum_hourly",Int16,
               ("longitude","latitude","hour"),attrib=att_var[1]);
    v.var[:] = data[3][:,:,1:nt,end];

    v = defVar(ds,"domain_yearly_minimum_hourly",Int16,
               ("longitude","latitude","hour"),attrib=att_var[1]);
    v.var[:] = data[4][:,:,1:nt,end];

    v = defVar(ds,"domain_yearly_mean_diurnalvariance",Int16,
           ("longitude","latitude"),attrib=att_var[2]);
    v.var[:] = data[1][:,:,nt+2,end];

    v = defVar(ds,"domain_yearly_std_diurnalvariance",Int16,
               ("longitude","latitude"),attrib=att_var[2]);
    v.var[:] = data[2][:,:,nt+2,end];

    v = defVar(ds,"domain_yearly_maximum_diurnalvariance",Int16,
               ("longitude","latitude"),attrib=att_var[2]);
    v.var[:] = data[3][:,:,nt+2,end];

    v = defVar(ds,"domain_yearly_minimum_diurnalvariance",Int16,
           ("longitude","latitude"),attrib=att_var[2]);
    v.var[:] = data[4][:,:,nt+2,end];

    v = defVar(ds,"domain_monthly_mean_climatology",Int16,
               ("longitude","latitude","month"),attrib=att_var[1]);
    v.var[:] = data[1][:,:,nt+1,1:12];

    v = defVar(ds,"domain_monthly_std_climatology",Int16,
               ("longitude","latitude","month"),attrib=att_var[1]);
    v.var[:] = data[2][:,:,nt+1,1:12];

    v = defVar(ds,"domain_monthly_maximum_climatology",Int16,
               ("longitude","latitude","month"),attrib=att_var[1]);
    v.var[:] = data[3][:,:,nt+1,1:12];

    v = defVar(ds,"domain_monthly_minimum_climatology",Int16,
               ("longitude","latitude","month"),attrib=att_var[1]);
    v.var[:] = data[4][:,:,nt+1,1:12];

    v = defVar(ds,"domain_monthly_mean_hourly",Int16,
               ("longitude","latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = data[1][:,:,1:nt,1:12];

    v = defVar(ds,"domain_monthly_std_hourly",Int16,
               ("longitude","latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = data[2][:,:,1:nt,1:12];

    v = defVar(ds,"domain_monthly_maximum_hourly",Int16,
               ("longitude","latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = data[3][:,:,1:nt,1:12];

    v = defVar(ds,"domain_monthly_minimum_hourly",Int16,
               ("longitude","latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = data[4][:,:,1:nt,1:12];

    v = defVar(ds,"domain_monthly_mean_diurnalvariance",Int16,
               ("longitude","latitude","month"),attrib=att_var[2]);
    v.var[:] = data[1][:,:,nt+2,1:12];

    v = defVar(ds,"domain_monthly_std_diurnalvariance",Int16,
               ("longitude","latitude","month"),attrib=att_var[2]);
    v.var[:] = data[2][:,:,nt+2,1:12];

    v = defVar(ds,"domain_monthly_maximum_diurnalvariance",Int16,
               ("longitude","latitude","month"),attrib=att_var[2]);
    v.var[:] = data[3][:,:,nt+2,1:12];

    v = defVar(ds,"domain_monthly_minimum_diurnalvariance",Int16,
               ("longitude","latitude","month"),attrib=att_var[2]);
    v.var[:] = data[4][:,:,nt+2,1:12];

    v = defVar(ds,"zonalavg_yearly_mean_climatology",Int16,
               ("latitude",),attrib=att_var[1]);
    v.var[:] = zdata[1][:,nt+1,end];

    v = defVar(ds,"zonalavg_yearly_std_climatology",Int16,
               ("latitude",),attrib=att_var[1]);
    v.var[:] = zdata[2][:,nt+1,end];

    v = defVar(ds,"zonalavg_yearly_maximum_climatology",Int16,
               ("latitude",),attrib=att_var[1]);
    v.var[:] = zdata[3][:,nt+1,end];

    v = defVar(ds,"zonalavg_yearly_minimum_climatology",Int16,
               ("latitude",),attrib=att_var[1]);
    v.var[:] = zdata[4][:,nt+1,end];

    v = defVar(ds,"zonalavg_yearly_mean_hourly",Int16,
               ("latitude","hour"),attrib=att_var[1]);
    v.var[:] = zdata[1][:,1:nt,end];

    v = defVar(ds,"zonalavg_yearly_std_hourly",Int16,
               ("latitude","hour"),attrib=att_var[1]);
    v.var[:] = zdata[2][:,1:nt,end];

    v = defVar(ds,"zonalavg_yearly_maximum_hourly",Int16,
               ("latitude","hour"),attrib=att_var[1]);
    v.var[:] = zdata[3][:,1:nt,end];

    v = defVar(ds,"zonalavg_yearly_minimum_hourly",Int16,
               ("latitude","hour"),attrib=att_var[1]);
    v.var[:] = zdata[4][:,1:nt,end];

    v = defVar(ds,"zonalavg_yearly_mean_diurnalvariance",Int16,
               ("latitude",),attrib=att_var[2]);
    v.var[:] = zdata[1][:,nt+2,end];

    v = defVar(ds,"zonalavg_yearly_std_diurnalvariance",Int16,
               ("latitude",),attrib=att_var[2]);
    v.var[:] = zdata[2][:,nt+2,end];

    v = defVar(ds,"zonalavg_yearly_maximum_diurnalvariance",Int16,
               ("latitude",),attrib=att_var[2]);
    v.var[:] = zdata[3][:,nt+2,end];

    v = defVar(ds,"zonalavg_yearly_minimum_diurnalvariance",Int16,
               ("latitude",),attrib=att_var[2]);
    v.var[:] = zdata[4][:,nt+2,end];

    v = defVar(ds,"zonalavg_monthly_mean_climatology",Int16,
               ("latitude","month"),attrib=att_var[1]);
    v.var[:] = zdata[1][:,nt+1,1:12];

    v = defVar(ds,"zonalavg_monthly_std_climatology",Int16,
               ("latitude","month"),attrib=att_var[1]);
    v.var[:] = zdata[2][:,nt+1,1:12];

    v = defVar(ds,"zonalavg_monthly_maximum_climatology",Int16,
               ("latitude","month"),attrib=att_var[1]);
    v.var[:] = zdata[3][:,nt+1,1:12];

    v = defVar(ds,"zonalavg_monthly_minimum_climatology",Int16,
               ("latitude","month"),attrib=att_var[1]);
    v.var[:] = zdata[4][:,nt+1,1:12];

    v = defVar(ds,"zonalavg_monthly_mean_hourly",Int16,
               ("latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = zdata[1][:,1:nt,1:12];

    v = defVar(ds,"zonalavg_monthly_std_hourly",Int16,
               ("latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = zdata[2][:,1:nt,1:12];

    v = defVar(ds,"zonalavg_monthly_maximum_hourly",Int16,
               ("latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = zdata[3][:,1:nt,1:12];

    v = defVar(ds,"zonalavg_monthly_minimum_hourly",Int16,
               ("latitude","hour","month"),attrib=att_var[1]);
    v.var[:] = zdata[4][:,1:nt,1:12];

    v = defVar(ds,"zonalavg_monthly_mean_diurnalvariance",Int16,
               ("latitude","month"),attrib=att_var[2]);
    v.var[:] = zdata[1][:,nt+2,1:12];

    v = defVar(ds,"zonalavg_monthly_std_diurnalvariance",Int16,
               ("latitude","month"),attrib=att_var[2]);
    v.var[:] = zdata[2][:,nt+2,1:12];

    v = defVar(ds,"zonalavg_monthly_maximum_diurnalvariance",Int16,
               ("latitude","month"),attrib=att_var[2]);
    v.var[:] = zdata[3][:,nt+2,1:12];

    v = defVar(ds,"zonalavg_monthly_minimum_diurnalvariance",Int16,
               ("latitude","month"),attrib=att_var[2]);
    v.var[:] = zdata[4][:,nt+2,1:12];

    v = defVar(ds,"meridionalavg_yearly_mean_climatology",Int16,
               ("longitude",),attrib=att_var[1]);
    v.var[:] = mdata[1][:,nt+1,end];

    v = defVar(ds,"meridionalavg_yearly_std_climatology",Int16,
               ("longitude",),attrib=att_var[1]);
    v.var[:] = mdata[2][:,nt+1,end];

    v = defVar(ds,"meridionalavg_yearly_maximum_climatology",Int16,
               ("longitude",),attrib=att_var[1]);
    v.var[:] = mdata[3][:,nt+1,end];

    v = defVar(ds,"meridionalavg_yearly_minimum_climatology",Int16,
               ("longitude",),attrib=att_var[1]);
    v.var[:] = mdata[4][:,nt+1,end];

    v = defVar(ds,"meridionalavg_yearly_mean_hourly",Int16,
               ("longitude","hour"),attrib=att_var[1]);
    v.var[:] = mdata[1][:,1:nt,end];

    v = defVar(ds,"meridionalavg_yearly_std_hourly",Int16,
               ("longitude","hour"),attrib=att_var[1]);
    v.var[:] = mdata[2][:,1:nt,end];

    v = defVar(ds,"meridionalavg_yearly_maximum_hourly",Int16,
               ("longitude","hour"),attrib=att_var[1]);
    v.var[:] = mdata[3][:,1:nt,end];

    v = defVar(ds,"meridionalavg_yearly_minimum_hourly",Int16,
               ("longitude","hour"),attrib=att_var[1]);
    v.var[:] = mdata[4][:,1:nt,end];

    v = defVar(ds,"meridionalavg_yearly_mean_diurnalvariance",Int16,
               ("longitude",),attrib=att_var[2]);
    v.var[:] = mdata[1][:,nt+2,end];

    v = defVar(ds,"meridionalavg_yearly_std_diurnalvariance",Int16,
               ("longitude",),attrib=att_var[2]);
    v.var[:] = mdata[2][:,nt+2,end];

    v = defVar(ds,"meridionalavg_yearly_maximum_diurnalvariance",Int16,
               ("longitude",),attrib=att_var[2]);
    v.var[:] = mdata[3][:,nt+2,end];

    v = defVar(ds,"meridionalavg_yearly_minimum_diurnalvariance",Int16,
               ("longitude",),attrib=att_var[2]);
    v.var[:] = mdata[4][:,nt+2,end];

    v = defVar(ds,"meridionalavg_monthly_mean_climatology",Int16,
               ("longitude","month"),attrib=att_var[1]);
    v.var[:] = mdata[1][:,nt+1,1:12];

    v = defVar(ds,"meridionalavg_monthly_std_climatology",Int16,
               ("longitude","month"),attrib=att_var[1]);
    v.var[:] = mdata[2][:,nt+1,1:12];

    v = defVar(ds,"meridionalavg_monthly_maximum_climatology",Int16,
               ("longitude","month"),attrib=att_var[1]);
    v.var[:] = mdata[3][:,nt+1,1:12];

    v = defVar(ds,"meridionalavg_monthly_minimum_climatology",Int16,
               ("longitude","month"),attrib=att_var[1]);
    v.var[:] = mdata[4][:,nt+1,1:12];

    v = defVar(ds,"meridionalavg_monthly_mean_hourly",Int16,
               ("longitude","hour","month"),attrib=att_var[1]);
    v.var[:] = mdata[1][:,1:nt,1:12];

    v = defVar(ds,"meridionalavg_monthly_std_hourly",Int16,
               ("longitude","hour","month"),attrib=att_var[1]);
    v.var[:] = mdata[2][:,1:nt,1:12];

    v = defVar(ds,"meridionalavg_monthly_maximum_hourly",Int16,
               ("longitude","hour","month"),attrib=att_var[1]);
    v.var[:] = mdata[3][:,1:nt,1:12];

    v = defVar(ds,"meridionalavg_monthly_minimum_hourly",Int16,
               ("longitude","hour","month"),attrib=att_var[1]);
    v.var[:] = mdata[4][:,1:nt,1:12];

    v = defVar(ds,"meridionalavg_monthly_mean_diurnalvariance",Int16,
               ("longitude","month"),attrib=att_var[2]);
    v.var[:] = mdata[1][:,nt+2,1:12];

    v = defVar(ds,"meridionalavg_monthly_std_diurnalvariance",Int16,
               ("longitude","month"),attrib=att_var[2]);
    v.var[:] = mdata[2][:,nt+2,1:12];

    v = defVar(ds,"meridionalavg_monthly_maximum_diurnalvariance",Int16,
               ("longitude","month"),attrib=att_var[2]);
    v.var[:] = mdata[3][:,nt+2,1:12];

    v = defVar(ds,"meridionalavg_monthly_minimum_diurnalvariance",Int16,
               ("longitude","month"),attrib=att_var[2]);
    v.var[:] = mdata[4][:,nt+2,1:12];

    close(ds);

    @info "$(Dates.now()) - Analysed $(info["source"]) $(info["product"]) data for the year $yr in $(gregionfullname(region)) has been saved into file $(fnc) and moved to the data directory $(fol)."

end


text:
module TestVine

using Copulas
using Base.Test

## equSets
##--------

@test Copulas.equSets([1, 2, 3, 4], [2, 1, 4, 3])
@test Copulas.equSets([1, 2, 2], [2, 1])

## pairMatch
##----------

## does match exist between two trees?

tr1 = Copulas.Tree(3, [1, 2], [4, 5, 6])
tr2 = Copulas.Tree(7, [2, 1], [4], [5, 6])
@test Copulas.pairMatch(tr1, tr2, 3) == [1 1]

tr2 = Copulas.Tree(7, [1], [1, 2], [5, 6])
@test Copulas.pairMatch(tr1, tr2, 3) == [1 2]

tr1 = Copulas.Tree(3, [1, 2], [4, 5])
tr2 = Copulas.Tree(7, [1], [1, 2], [5, 4])
@test Copulas.pairMatch(tr1, tr2, 3) == [1 2; 2 3]

tr1 = Copulas.Tree(3, [1, 2])
tr2 = Copulas.Tree(7, [3, 4])
@test isempty(Copulas.pairMatch(tr1, tr2, 3))

## findPairMatch
##--------------

## for array of trees, where do matches exist?
## also: incomplTrees tests

trs = [Copulas.Tree(1, [2], [3]),
       Copulas.Tree(2, [1, 3]),
       Copulas.Tree(3, [1, 2])
       ]

@test isempty(Copulas.incomplTrees(trs))

trs = [Copulas.Tree(1, [3]);
       Copulas.Tree(2, [3], [5]);
       Copulas.Tree(3, [1], [2]);
       Copulas.Tree(4, [5], [6]);
       Copulas.Tree(5, [2], [4]);
       Copulas.Tree(6, [4])]

@test Copulas.incomplTrees(trs) == [1, 2, 3, 4, 5, 6]

expMatches = [Copulas.CondEdge(1, 2, [3], [3]),
              Copulas.CondEdge(2, 4, [5], [5]),
              Copulas.CondEdge(3, 5, [2], [2]),
              Copulas.CondEdge(5, 6, [4], [4])]

@test expMatches == Copulas.findPairMatch(trs, 2, Copulas.incomplTrees(trs))

## vertexTable
##------------

expVTable = Array{Int, 1}[[1],
                          [1, 2],
                          [3],
                          [2],
                          [3, 4],
                          [4]]
@test expVTable == Copulas.vertexTable(expMatches, 6)

## findPairMatch
##--------------

## given first dvine layer
trs = [Copulas.Tree(1, [2]),
       Copulas.Tree(2, [1], [3]),
       Copulas.Tree(3, [2], [4]),
       Copulas.Tree(4, [3])]

expMatches = [Copulas.CondEdge(1, 3, [2], [2]),
              Copulas.CondEdge(2, 4, [3], [3])]

foundPairs = Copulas.findPairMatch(trs, 2, Copulas.incomplTrees(trs))
@test expMatches == foundPairs

## reachedVertices
##----------------

@test [1, 2, 3, 4] == Copulas.reachedVertices(foundPairs)

edges = [Copulas.CondEdge(1, 3, [2], [2]),
         Copulas.CondEdge(1, 4, [3], [3])]

@test [1, 3, 4] == Copulas.reachedVertices(edges)

## mandatoryEdges
##---------------

edges = [Copulas.CondEdge(1, 3, [2], [2]),
         Copulas.CondEdge(2, 4, [3], [3])]
verticesToBeLinked = [1, 2, 3, 4]

expChosen = Bool[1, 1]
actChosen = Copulas.mandatoryEdges(edges, verticesToBeLinked, 4)[1]
@test expChosen == actChosen

expRemainLinks = []
actRemainLinks = Copulas.mandatoryEdges(edges, verticesToBeLinked, 4)[2]
@test expRemainLinks == actRemainLinks

edges = [Copulas.CondEdge(1, 3, [2], [2]),
         Copulas.CondEdge(2, 3, [4], [4]),
         Copulas.CondEdge(2, 4, [3], [3]),
         Copulas.CondEdge(3, 4, [2], [2])]
verticesToBeLinked = [1, 2, 3, 4]

expChosen = Bool[1, 0, 0, 0]
actChosen = Copulas.mandatoryEdges(edges, verticesToBeLinked, 4)[1]
@test expChosen == actChosen

expRemainLinks = [2, 4]
actRemainLinks = Copulas.mandatoryEdges(edges, verticesToBeLinked, 4)[2]
@test expRemainLinks == actRemainLinks

## selectLinks
##------------

## second layer d-Vine
trs = [Copulas.Tree(1, [2]),
       Copulas.Tree(2, [1], [3]),
       Copulas.Tree(3, [2], [4]),
       Copulas.Tree(4, [3])]

actLinks = Copulas.selectLinks(trs, 2)
expLinks = [Copulas.CondEdge(1, 3, [2], [2]),
            Copulas.CondEdge(2, 4, [3], [3])]
@test actLinks == expLinks

## third layer d-Vine
##-------------------

trs = [Copulas.Tree(1, [2, 3]),
       Copulas.Tree(2, [1], [3, 4]),
       Copulas.Tree(3, [2, 1], [4]),
       Copulas.Tree(4, [3, 2])]

actLinks = Copulas.selectLinks(trs, 3)
expLinks = [Copulas.CondEdge(1, 4, [2, 3], [3, 2])]
@test actLinks == expLinks

## test DVine construction
##------------------------

## expected result
dVn4 = Copulas.Vine([0 2 2 2;
                   1 0 3 3;
                   2 2 0 4;
                   3 3 3 0])

## given first dvine layer
trs = [Copulas.CTreePaths(1, [2]),
       Copulas.CTreePaths(2, [1], [3]),
       Copulas.CTreePaths(3, [2], [4]),
       Copulas.CTreePaths(4, [3])]
vn = Copulas.Vine(trs)

## construct D-Vine
vnConstructed = Copulas.vineFrom2(vn)

@test vnConstructed == dVn4

end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 2d679c48-8af5-11eb-32fc-c9a76a667e61
begin
	using OnlineStats, Plots, Random
	theme(:lime)
end

# ╔═╡ 7ac57a82-8aff-11eb-0bbb-9fc037be9621
md"![](https://user-images.githubusercontent.com/8075494/111925031-87462b80-8a7d-11eb-98e2-eae044b13a3f.png)"

# ╔═╡ 497276a6-8af5-11eb-0d24-3d6c0c01c523
n = 100

# ╔═╡ 5ac71308-8af5-11eb-3767-2b7ef25271e1
nframes = 100

# ╔═╡ efb59472-8b03-11eb-17b7-d3c7960eb49f
begin
	o = HeatMap(-5:.2:5, 0:.2:10)

	anim = @animate for i in 1:nframes
		x = randn(5i)
		y = randexp(5i)
		fit!(o, zip(x,y))
		plot(o)
	end
	gif(anim, "temp.gif", fps=10)
end

# ╔═╡ Cell order:
# ╟─7ac57a82-8aff-11eb-0bbb-9fc037be9621
# ╟─2d679c48-8af5-11eb-32fc-c9a76a667e61
# ╠═497276a6-8af5-11eb-0d24-3d6c0c01c523
# ╠═5ac71308-8af5-11eb-3767-2b7ef25271e1
# ╠═efb59472-8b03-11eb-17b7-d3c7960eb49f


text:
import Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using Revise

using Beamforming
using Plots
import Seis, Seis.Plot
using BenchmarkTools
import FFTW

function shift!(s::Seis.Trace, Δ)
    s.t .= circshift(s.t, round(Int, Δ/s.delta))
    s
end

shift!(s::AbstractArray{<:Seis.Trace}, Δ) = shift!.(s, Δ)

"Cyclically shift a SAC trace backward by Δ s in time"
function fdshift!(s::Seis.Trace{T}, Δ) where T
    S = FFTW.fft(Seis.trace(s))
    δ = s.delta
    nf = Int32(length(S))
    coeff = 2π*im*Δ/(δ * 2*(nf - 1))
    @inbounds for i in 1:nf
        S[i] = S[i]*exp(coeff*(i-1))
    end
    s.t .= real.(FFTW.ifft(S))
    s
end

function fdshift!(s::AbstractArray{<:Seis.Trace}, Δ::AbstractArray)
    @assert length(s) == length(Δ)
    S = FFTW.fft(hcat(Seis.trace.(s)...), 1)
    δ = s[1].delta
    nf, nsta = size(S)
    coeff = Complex{Float32}(2π*im/(δ*2*(nf-1)))
    @inbounds for j in 1:nsta, i in 1:nf
        S[i,j] = S[i,j]*exp(-coeff*Float32(Δ[j])*(i - 1))
    end
    FFTW.ifft!(S, 1)
    @inbounds for j in 1:nsta
        s[j].t .= real.(S[:,j])
    end
    s
end

# using ArrayFire
# function fdshift_af!(s::SACtr, Δ)
#     δ = s[:delta]
#     S = AFArray(Complex{Float32}.(s.t))
#     fft!(S)
#     nf = length(S)
#     coeffs = AFArray(Complex{Float32}.(2π.*im.*Δ./δ.*collect(0:(nf-1))./(2*(nf-1))))
#     S .= S.*exp(coeffs)
#     ifft!(S)
#     s[:t] = real(S)
#     s
# end
#
# function fdshift_af!(s::AbstractArray{SACtr}, Δ::AbstractArray)
#     @assert length(s) == length(Δ)
#     @assert all(s[:npts] .== s[1][:npts])
#     @assert all(isapprox.(s[:b], s[1][:b], atol=s[1][:delta]/4))
#     δ = s[1][:delta]
#     nsta = length(s)
#     npts = s[1][:npts]
#     np = nextpow2(s[1][:npts])
#     nf = np
#     S = zeros(Complex{Float32}, np, nsta)
#     S[1:npts,:] .= hcat(s[:t]...)
#     S = AFArray(S)
#     fft!(S)
#     coeffs = AFArray(Complex{Float32}.(2π.*im.*Δ'./δ.*collect(0:(nf-1))./(2*(nf-1))))
#     S = S.*exp(coeffs)
#     S = Array(real(ifft!(S, 1)))
#     @inbounds for i in 1:nsta
#         s[i].t[:] .= S[1:npts,i]
#     end
#     SAC.update_headers!(s)
# end

function beamformfd(s::AbstractArray{<:Seis.Trace}, t1, t2, sx1, sx2, sy1, sy2, ds)
    nsta = length(s)
    s = Seis.cut.(s, t1, t2)
    npts = Seis.nsamples(s[1])
    @assert all(Seis.nsamples.(s) .== npts)
    δ = s[1].delta
    sxs, sys = sx1:ds:sx2, sy1:ds:sy2
    sx = Beamforming.s_per_km.(sxs)
    sy = Beamforming.s_per_km.(sys)
    slow = hcat(vec(sxs), vec(sys))
    nslow = size(slow, 1)
    pow = zeros(Float32, nslow)
    Δ = Beamforming.delay_times(s, sx, sy)
    spec = FFTW.fft(hcat(Seis.trace.(s)...), 1)
    nf = size(spec, 1)
    # Steering vectors
    steers = 2π*im*collect(0:(nf-1))./(δ*2*(nf - 1))
    steer = Array{Complex{Float32}}(undef, nf, nsta)
    shift_spec = similar(spec)
    for k in eachindex(sys), j in eachindex(sxs)
        for i in 1:nf
            shift_spec[i,:] .= spec[i,:].*exp.(steers[i].*Δ[:,j,k])
        end
        for ista in 1:nsta
            temp = zero(eltype(shift_spec))
            for ifreq in 1:nf
                temp += shift_spec[ifreq,ista]
            end
            pow[j,k] += abs(temp^2)
        end
    end
    sxs, sys, pow
end

# function beamform_af(s::AbstractArray{SACtr}, Δ::AbstractArray, slow)
#     @assert length(s) == length(Δ)
#     @assert all(s[:npts] .== s[1][:npts])
#     @assert all(isapprox.(s[:b], s[1][:b], atol=s[1][:delta]/4))
#     δ = s[1][:delta]
#     nsta = length(s)
#     npts = s[1][:npts]
#     np = nextpow2(s[1][:npts])
#     nf = np
#     S = zeros(Complex{Float32}, np, nsta)
#     S[1:npts,:] .= hcat(s[:t]...)
#     S = AFArray(S)
# end

s = Seis.sample_data(:array)
Seis.cut!.(s, maximum(s.b), minimum(Seis.endtime.(s)))
Seis.cut!.(s, s.b, s.b.+s[1].delta.*(Seis.nsamples(s[1])-1))
# @btime fdshift!($s, $(-s[:a]))
# @btime fdshift_af!($s, $(-s[:a]))
@time s′ = shift!(deepcopy(s), -s.picks.A.time)
@time s″ = fdshift!(deepcopy(s), -2s.picks.A.time)

@time beamformfd(s, 1050, 1200, -3, 3, -3, 3, 0.5)
@time Beamforming.beamform(s, 1050, 1200, -3, 3, -3, 3, 0.5)

# Plots.plot([Seis.Plot.section(ss, align=ss.b) for ss in (s, s′, s″)]...)

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
abstract Report

function balancetable(accounts)

end

immutable BSSection
    heading::Compat.UTF8String
    subsections::Vector{BSSection}
    entries::Vector{Pair{Compat.UTF8String, Basket}}

    BSSection(heading) = new(heading, [], [])
end

function total(s::BSSection)
    acc = zero(Basket)
    for sect in s.subsections
        acc += total(sect)
    end
    for (_, account) in s.entries
        acc += account
    end
    acc
end

function Base.push!{T<:AbstractString}(s::BSSection, p::Pair{T, Basket})
    push!(s.entries, Compat.UTF8String(p[1]) => p[2])
end

function Base.push!(s::BSSection, p::BSSection)
    push!(s.subsections, p)
end

immutable BalanceSheet <: Report
    assets::BSSection
    le::BSSection
end

function balancesheet(ledger::Ledger)
    bals = balances(ledger)
    assets = BSSection("Assets")
    liabilities = BSSection("Liabilities")
    equity = BSSection("Equity")
    le = BSSection("Liabilities & Equity")
    earnings = Basket()
    unrealized = Basket()
    for (acc, bal) in bals
        @match accounttype(acc) begin
            :asset => push!(assets, name(acc) => bal)
            :liability => push!(liabilities, name(acc) => -bal)
            :equity => push!(equity, name(acc) => -bal)
            :revenue => push!(earnings, -bal)
            :expense => push!(earnings, -bal)
            :trading => push!(unrealized, -bal)
        end
    end
    push!(le, liabilities)
    push!(le, equity)
    push!(equity, "Retained Earnings" => earnings)
    push!(equity, "Unrealized Gains" => unrealized)
    BalanceSheet(assets, le)
end

Currencies.valuate(table, as::Symbol, p::Pair{String, Basket}) =
    p.first => Basket(valuate(table, as, p.second))

function Currencies.valuate(table, as::Symbol, s::BSSection)
    ns = BSSection(s.heading)
    for ss in s.subsections
        push!(ns, valuate(table, as, ss))
    end
    for ent in s.entries
        push!(ns, valuate(table, as, ent))
    end
    ns
end

Currencies.valuate(table, as::Symbol, bs::BalanceSheet) =
    BalanceSheet(valuate(table, as, bs.assets), valuate(table, as, bs.le))


text:
function driveHex(fgl, posecount::Int; steps::Int=5)
    # Drive around in a hexagon
    for i in (posecount-1):(posecount-1+steps)
        psym = Symbol("x$i")
        posecount += 1
        nsym = Symbol("x$(i+1)")
        addVariable!(fgl, nsym, Pose2)
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
        addFactor!(fgl, [psym;nsym], pp, autoinit=false )
    end

    return posecount
end

function offsetHexLeg(dfg::G, posecount::Int; direction=:right) where G <: AbstractDFG
    psym = Symbol("x$(posecount-1)")
    nsym = Symbol("x$(posecount)")
    posecount += 1
    addVariable!(dfg, nsym, Pose2)
    pp = nothing
    if direction == :right
        pp = Pose2Pose2(MvNormal([10.0;0;-pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    elseif direction == :left
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    end
    addFactor!(dfg, [psym; nsym], pp, autoinit=false )
    return posecount
end


text:
#---test emulating for loop , i, upperbound, lenA +s - 1
a = [2, 3, 4]
total = length(a)
#correct for s = 1
#------for 0
s = 1
i = s #initial i0

#-----start loop
i # <----
#-------- variables 3PICK WHICH ONE Equals s
#_final = copy(i - s + 1) #_len - s + i # ERROR: 4 #issue i + s #the Most  correct!
_start1 = copy(i + s) #2 (should be 1) #4 (array out of bounds error )
_start2 = copy(i + s - 1) #working well #as a start  #1  #last = 3  #inbounds #accepted

_final6 = copy(length(a) + s - 1) # works well, for a finish #3 #inbounds #accepted
_final6_5 = copy(length(a) + s + i - 1)
_final7 = copy(length(a) + s + i - 3) #4
_final7 = copy(s + i) # only correct   #4

_final8 = copy(s + i - 1) # works well as an end  #3 #least constraints #adopt this
#infer formula: s + i -1 if i = 1 : s +1 -1 = s (start init value ): either 0 (0th index) or 1 (1st index)
#-----------------------# seems you plugged in i = s - 1
i += 1
#-----------------------
lenA = length(a)
lenA + i
lenA + s - 1 + -lenA - 1  # here , if s = i = 1
lenA + s - 1 + (lenA + 1) - lenA - 1 - (lenA + s - 1) #equivalent (already) #or by composition
lenA(1 - 1) + s - 2  #experimental # same #equivalent
s - 2 # here, is start at s=i=1 it yeilds -1 # not acceptable in the positive domain
# mind boggling line abve seems to have no equal (initially )
#-----------adopt this for for every future endevour
#s either 0 or 1
for i ∈ 1:lenA+s-1  # idea= keep fixating starting i at  i=1 #strict condition: needs s=i=1 to work well
    println(i) #prints range correctly 1,2,3 (accordingly)
end

#testing ranges
s = 1
i = 1 #initial i0 #i is fixed to 1

#-----start loop
i # <----
#-------- variables 3PICK WHICH ONE Equals s
#_final = copy(i - s + 1) #_len - s + i # ERROR: 4 #issue i + s #the Most  correct!
_start1 = copy(i + s)
_start2 = copy(i + s - 1)
#_final = copy(i + s - 1) #_len - s + i # ERROR: 4 #issue i + s #(al)most correct
#_final4 = copy(total - s + 1) # i believe so
#_final5 = copy(length(a) - 1)
_final6 = copy(length(a) + s - 1)
_final7 = copy(length(a) + s + i - 3)
_final7 = copy(s + i) # only correct

_final8 = copy(s + i - 1)
#-----------------------
i += 1
#------- fixed values [Desired]
_first = s #s
_len = calcLength(2, 4; s) # 2, 3, 4 #euclidDist + start= 1  # |b - a| + s # always 3
lenA = copy(length(a)) # always 3
_endingIndex = _len   # + i - s #value inside array bracket
_len - lenA #always 0
_len === lenA

_finalAlt = copy(i + s) #seems ok  #errorneous
_final2 = copy(total - s) # i believe so
_final3 = copy(total - s + i) # i believe so
#----------------

#metric:
_endingIndex == lenA
_endingIndex - lenA + i #start

a[_len]
#recall i=0 ; length - 1
#when   i=1 ; length -1

# endingIndex == _len - s

#infer

#best start

#best final

_final = lenA +s +i -1

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
include("context.jl")
include("plugins.jl")
include("superchunk.jl")
include("metalayers.jl")


text:
# store psis in vector and iterate to do derivatives

function solovev_psi_functions()
    # up down symmetric homogeneous solutions
    psi1(x,y) = 1
    psi2(x,y) = x^2
    psi3(x,y) = y^2 - x^2*log(x)
    psi4(x,y) = x^4 - 4*x^2*y^2
    psi5(x,y) = 2*y^4 - 9*y^2*x^2 + 3*x^4*log(x) - 12*x^2*y^2*log(x)
    psi6(x,y) = x^6 - 12*x^4*y^2 + 8*x^2*y^4
    psi7(x,y) = 8*y^6 - 140*y^4*x^2 + 75*y^2*x^4 - 15*x^6*log(x) + 180*x^4*y^2*log(x) - 120*x^2*y^4*log(x)

    # up down asymmetric homogenous solutions
    psi8(x,y) = y
    psi9(x,y) = y*x^2
    psi10(x,y) = y^3 - 3*y*x^2*log(x)
    psi11(x,y) = 3*y*x^4 - 4*y^3*x^2
    psi12(x,y) = 8*y^5 - 45*y*x^4 - 80*y^3*x^2*log(x) + 60*y*x^4*log(x)

    psibasis = [psi1 psi2 psi3 psi4 psi5 psi6 psi7 psi8 psi9 psi10 psi11 psi12]

    return psibasis
end

function solve_solovev_psi_functions(x, y) #XXX might not want to call this function the same as the function above?
    # up down symmetric homogeneous solutions
    psi1 = 1
    psi2 = x^2
    psi3 = y^2 - x^2*log(x)
    psi4 = x^4 - 4*x^2*y^2
    psi5 = 2*y^4 - 9*y^2*x^2 + 3*x^4*log(x) - 12*x^2*y^2*log(x)
    psi6 = x^6 - 12*x^4*y^2 + 8*x^2*y^4
    psi7 = 8*y^6 - 140*y^4*x^2 + 75*y^2*x^4 - 15*x^6*log(x) + 180*x^4*y^2*log(x) - 120*x^2*y^4*log(x)

    # up down asymmetric homogenous solutions
    psi8 = y
    psi9 = y*x^2
    psi10 = y^3 - 3*y*x^2*log(x)
    psi11 = 3*y*x^4 - 4*y^3*x^2
    psi12 = 8*y^5 - 45*y*x^4 - 80*y^3*x^2*log(x) + 60*y*x^4*log(x)

    psivals = [psi1 psi2 psi3 psi4 psi5 psi6 psi7 psi8 psi9 psi10 psi11 psi12]

    return psivals
end

function solovev_parameterize(n, delta, epsilon, kappa)
    # epsilon = reverse aspect ratio (epsilon = a/R0 for R0 the major radius (center of machine to
    # center of region of plasma) and a the minor radius (center of plasma to wall))
    # kappa = elongation
    # sin(alpha) = delta = triangularity
    # n = resolution of points in x/y

    alpha = asin(delta)

    # calculating x and y coordinates of flux surface cross section
    tau = LinRange(0, 2*pi, n)
    x = zeros(n)
    y = zeros(n)
    for i=1:n
        x[i] = 1 + epsilon*cos(tau[i] + alpha*sin(tau[i]))
        y[i] = epsilon*kappa*sin(tau[i])
    end

    return tau, x, y
end

function solovev_particular_psi(A, x, y)
    psi = x^4/8 + A * (1/2*x^2*log(x) - x^4/8)
    return psi
end

function solovev_homo_partial_x(functs, xval, yval) #XXX could be confusing that the input variable psis is also called the same as two seperate functions defined above.
    partial_x = zeros(12)
    for i=1:12
        # partial_x[i] = ForwardDiff.derivative(x -> psis[i](x,y), x)

        # XXX: ForwardDiff.derivative returns a number, not a function. The input looks like (::Function, ::Real), and the function has to be single variable
        #      so it doesn't know what to do with psis[i](x,y), it has to be psis[i](x,yval)
        partial_x[i] = ForwardDiff.derivative(x -> functs[i](x,yval), xval) #XXX <-- applying this correction to the other instances should fix the issue.
    end
    # for j=1:12
    #     partial_x[j] = partial_x[j][xval, yval]
    # end

    return partial_x
end

function solovev_homo_partial_xx(functs, xval, yval)
    # px = []
    pxx = zeros(12)
    for i=1:12
        # XXX Since ForwardDiff isn't returning functions, this is the way to take a second derivative:
        pxx[i] = ForwardDiff.derivative(x -> ForwardDiff.derivative(x -> functs[i](x,yval), x), xval)
    end
    return pxx
end

function solovev_homo_partial_y(functs, xval, yval) #XXX Need to just apply the changes above to the y versions.
    partial_y = zeros(12)
    for i=1:12
        partial_y[i] = ForwardDiff.derivative(y -> functs[i](xval,y), yval)
    end
    return partial_y
end

function solovev_homo_partial_yy(functs, xval, yval)
    # py = zeros(12)
    pyy = zeros(12)
    for i=1:12
        # py[i] = ForwardDiff.derivative(y -> psis[i](x,y), y)
        pyy[i] = ForwardDiff.derivative(y -> ForwardDiff.derivative(y -> functs[i](xval,y), y), yval)
    end
    return pyy
end

function solovev_particular_partial_x(A, xval, yval)
    partial_x = ForwardDiff.derivative(x -> solovev_particular_psi(A, x, yval), xval)
    return partial_x
end

function solovev_particular_partial_xx(A, xval, yval)
    # partial_x = ForwardDiff.derivative(x -> solovev_particular_psi(A, x, y), x)
    partial_xx = ForwardDiff.derivative(x -> ForwardDiff.derivative(x -> solovev_particular_psi(A, x, yval), x), xval)
    return partial_xx
end

function solovev_particular_partial_y(A, xval, yval)
    partial_y = ForwardDiff.derivative(y -> solovev_particular_psi(A, xval, y), yval)
    return partial_y
end

function solovev_particular_partial_yy(A, xval, yval)
    # partial_y = ForwardDiff.derivative(y -> solovev_particular_psi(A)(y), yval)
    partial_yy = ForwardDiff.derivative(y -> ForwardDiff.derivative(y -> solovev_particular_psi(A, xval, y), y), yval)
    return partial_yy
end

function solovev_flux_function(A, delta, epsilon, kappa, B0; upsep=[0,0], downsep=[0,0])
    alpha = asin(delta)

    # last terms are xsep/ysep - can be changed to be inputs if need be
    N1 = -(1 + alpha)^2/(epsilon*kappa^2)
    N2 = (1 - alpha)^2/(epsilon*kappa^2)
    N3 = -kappa/(epsilon*(cos(alpha)^2))

    x = [1 + epsilon, 1 - epsilon, 1 - delta*epsilon, upsep[1], downsep[1]]
    y = [0, 0, kappa*epsilon, upsep[2], downsep[2]]

    psisbasis = solovev_psi_functions()

    rhs1 = solovev_particular_psi(A, x[1], y[1])
    rhsy2 = solovev_particular_partial_y(A, x[1], y[1])
    rhsyy3 = solovev_particular_partial_yy(A, x[1], y[1])
    rhsx3 = solovev_particular_partial_x(A, x[1], y[1])

    rhs4 = solovev_particular_psi(A, x[2], y[2])
    rhsy5 = solovev_particular_partial_y(A, x[2], y[2])
    rhsyy6 = solovev_particular_partial_yy(A, x[2], y[2])
    rhsx6 = solovev_particular_partial_x(A, x[2], y[2])

    psis1 = solve_solovev_psi_functions(x[1], y[1])
    psisy2 = solovev_homo_partial_y(psisbasis, x[1], y[1])
    psisyy3 = solovev_homo_partial_yy(psisbasis, x[1], y[1])
    psisx3 = solovev_homo_partial_x(psisbasis, x[1], y[1])

    psis4 = solve_solovev_psi_functions(x[2], y[2])
    psisy5 = solovev_homo_partial_y(psisbasis, x[2], y[2])
    psisyy6 = solovev_homo_partial_yy(psisbasis, x[2], y[2])
    psisx6 = solovev_homo_partial_x(psisbasis, x[2], y[2])

    Amat = zeros(12, 12)
    Amat[1,:] = psis1
    Amat[2,:] = psisy2
    Amat[3,:] = psisyy3 + N1*psisx3
    Amat[4,:] = psis4
    Amat[5,:] = psisy5
    Amat[6,:] = psisyy6 + N2*psisx6

    if upsep==[0,0] && downsep==[0,0]
        # calculating values of flux function and derivatives at boundaries
        rhs7 = solovev_particular_psi(A, x[3], y[3])
        rhsx8 = solovev_particular_partial_x(A, x[3], y[3])
        rhsxx9 = solovev_particular_partial_xx(A, x[3], y[3])
        rhsy9 = solovev_particular_partial_y(A, x[3], y[3])
        rhs10 = solovev_particular_psi(A, x[3], -y[3])
        rhsx11 = solovev_particular_partial_x(A, x[3], -y[3])
        rhsxx12 = solovev_particular_partial_xx(A, x[3], -y[3])
        rhsy12 = solovev_particular_partial_y(A, x[3], -y[3])

        psis7 = solve_solovev_psi_functions(x[3], y[3])
        psisx8 = solovev_homo_partial_x(psisbasis, x[3], y[3])
        psisxx9 = solovev_homo_partial_xx(psisbasis, x[3], y[3])
        psisy9 = solovev_homo_partial_y(psisbasis, x[3], y[3])
        psis10 = solve_solovev_psi_functions(x[3], -y[3])
        psisx11 = solovev_homo_partial_x(psisbasis, x[3], -y[3])
        psisxx12 = solovev_homo_partial_xx(psisbasis, x[3], -y[3])
        psisy12 = solovev_homo_partial_y(psisbasis, x[3], -y[3])

        Amat[7,:] = psis7
        Amat[8,:] = psisx8
        Amat[9,:] = psisxx9 + N3*psisy9
        Amat[10,:] = psis10
        Amat[11,:] = psisx11
        Amat[12,:] = psisxx12 - N3*psisy12

        b = [rhs1, rhsy2, rhsyy3+N1*rhsx3, rhs4, rhsy5, rhsyy6+N2*rhsx6,
        rhs7, rhsx8, rhsxx9+N3*rhsy9, rhs10, rhsx11, rhsxx12-N3*rhsy12]
    elseif downsep==[0,0]
        rhs7 = solovev_particular_psi(A, x[4], y[4])
        rhsx8 = solovev_particular_partial_x(A, x[4], y[4])
        rhsy9 = solovev_particular_partial_y(A, x[4], y[4])
        rhs10 = solovev_particular_psi(A, x[3], -y[3])
        rhsx11 = solovev_particular_partial_x(A, x[3], -y[3])
        rhsxx12 = solovev_particular_partial_xx(A, x[3], -y[3])
        rhsy12 = solovev_particular_partial_y(A, x[3], -y[3])

        psis7 = solve_solovev_psi_functions(x[4], y[4])
        psisx8 = solovev_homo_partial_x(psisbasis, x[4], y[4])
        psisy9 = solovev_homo_partial_y(psisbasis, x[4], y[4])
        psis10 = solve_solovev_psi_functions(x[3], -y[3])
        psisx11 = solovev_homo_partial_x(psisbasis, x[3], -y[3])
        psisxx12 = solovev_homo_partial_xx(psisbasis, x[3], -y[3])
        psisy12 = solovev_homo_partial_y(psisbasis, x[3], -y[3])

        Amat[7,:] = psis7
        Amat[8,:] = psisx8
        Amat[9,:] = psisy9
        Amat[10,:] = psis10
        Amat[11,:] = psisx11
        Amat[12,:] = psisxx12 - N3*psisy12

        b = [rhs1, rhsy2, rhsyy3+N1*rhsx3, rhs4, rhsy5, rhsyy6+N2*rhsx6,
        rhs7, rhsx8, rhsy9, rhs10, rhsx11, rhsxx12+N3*rhsy12]
    elseif upsep==[0,0]
        rhs7 = solovev_particular_psi(A, x[3], y[3])
        rhsx8 = solovev_particular_partial_x(A, x[3], y[3])
        rhsxx9 = solovev_particular_partial_xx(A, x[3], y[3])
        rhsy9 = solovev_particular_partial_y(A, x[3], y[3])
        rhs10 = solovev_particular_psi(A, x[5], y[5])
        rhsx11 = solovev_particular_partial_x(A, x[5], y[5])
        rhsy12 = solovev_particular_partial_y(A, x[5], y[5])

        psis7 = solve_solovev_psi_functions(x[3], y[3])
        psisx8 = solovev_homo_partial_x(psisbasis, x[3], y[3])
        psisxx9 = solovev_homo_partial_xx(psisbasis, x[3], y[3])
        psisy9 = solovev_homo_partial_y(psisbasis, x[3], y[3])
        psis10 = solve_solovev_psi_functions(x[5], y[5])
        psisx11 = solovev_homo_partial_x(psisbasis, x[5], y[5])
        psisy12 = solovev_homo_partial_y(psisbasis, x[5], y[5])

        Amat[7,:] = psis7
        Amat[8,:] = psisx8
        Amat[9,:] = psisxx9 + N3*psisy9
        Amat[10,:] = psis10
        Amat[11,:] = psisx11
        Amat[12,:] = psisy12

        b = [rhs1, rhsy2, rhsyy3+N1*rhsx3, rhs4, rhsy5, rhsyy6+N2*rhsx6,
        rhs7, rhsx8, rhsxx9+N3*rhsy9, rhs10, rhsx11, rhsy12]
    else
        rhs7 = solovev_particular_psi(A, x[4], y[4])
        rhsx8 = solovev_particular_partial_x(A, x[4], y[4])
        rhsy9 = solovev_particular_partial_y(A, x[4], y[4])
        rhs10 = solovev_particular_psi(A, x[5], y[5])
        rhsx11 = solovev_particular_partial_x(A, x[5], y[5])
        rhsy12 = solovev_particular_partial_y(A, x[5], y[5])

        psis7 = solve_solovev_psi_functions(x[4], y[4])
        psisx8 = solovev_homo_partial_x(psisbasis, x[4], y[4])
        psisy9 = solovev_homo_partial_y(psisbasis, x[4], y[4])
        psis10 = solve_solovev_psi_functions(x[5], y[5])
        psisx11 = solovev_homo_partial_x(psisbasis, x[5], y[5])
        psisy12 = solovev_homo_partial_y(psisbasis, x[5], y[5])

        Amat[7,:] = psis7
        Amat[8,:] = psisx8
        Amat[9,:] = psisy9
        Amat[10,:] = psis10
        Amat[11,:] = psisx11
        Amat[12,:] = psisy12

        b = [rhs1, rhsy2, rhsyy3+N1*rhsx3, rhs4, rhsy5, rhsyy6+N2*rhsx6,
        rhs7, rhsx8, rhsy9, rhs10, rhsx11, rhsy12]
    end

    c = Amat \ -b
    psifunctions = solovev_psi_functions()
    psi(x,y) = solovev_particular_psi(A,x,y) .+ dot(c, [p(x,y) for p in psifunctions])

    Bx(x, y) = ForwardDiff.derivative(u -> -psi(x,u), y)
    By(x, y) = ForwardDiff.derivative(u -> psi(u,y), x) #XXX changed the sign of this to make consistent with trace.jl, should check later which way is actually correct.
    # B0 = R0^4*B0^2/psi0^2
    Bphi(x,y) = (B0^2 - 2*A*psi(x,y))^(1/2)
    B(x,y) = norm(Float64[Bx(x,y), By(x,y), Bphi(x,y)])

    bx(x,y) = Bx(x,y) / B(x,y)
    by(x,y) = By(x,y) / B(x,y)
    bz(x,y) = Bphi(x,y) / B(x,y)

    return psi, bx, by, bz
end


text:
using Documenter, MITgcmTools
import PlutoSliderServer
using Plots, CairoMakie

makedocs(;
    modules=[MITgcmTools],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Manual" => "functionalities.md",
        "Examples" => "examples.md",
#        "Examples" => Any[
#            "Guide " => "examples.md",
#            "Listing" => [map(s -> "generated/$(s[1:end-2])md",lst)...],
#            ],
    ],
    repo="https://github.com/gaelforget/MITgcmTools.jl/blob/{commit}{path}#L{line}",
    sitename="MITgcmTools.jl",
    authors="gaelforget <gforget@mit.edu>",
)

pth = joinpath(@__DIR__, "build","examples")
lst=("HS94_animation.jl","HS94_particles.jl","MITgcm_configurations.jl","MITgcm_run.jl","MITgcm_worklow.jl","HS94_Makie.jl","MITgcm_scan_output.jl")
for i in lst
    fil_in=joinpath(@__DIR__,"..","examples",i)
    fil_out=joinpath(pth,i[1:end-2]*"html")
    PlutoSliderServer.export_notebook(fil_in)
    mv(fil_in[1:end-2]*"html",fil_out)
    #cp(fil_in[1:end-2]*"html",fil_out)
end

deploydocs(;
    repo="github.com/gaelforget/MITgcmTools.jl",
)


text:
abstract type AbstractPooledDenseArray{T,N} <: DenseArray{T,N}
end

mutable struct PooledArray{T,N} <: AbstractPooledDenseArray{T,N}
    array::Array{T,N}
    ref_count::Int
    dispose::Function
    function PooledArray(a::Array{T,N}, dispose::Function) where {T,N}
        new{T,N}(a, 1, dispose)
    end
end

resource(a::AbstractPooledDenseArray{T,N}) where {T,N} = a.array

Base.size(a::AbstractPooledDenseArray{T,N}) where {T,N} = size(a.array)

Base.IndexStyle(::Type{<:AbstractPooledDenseArray{T,N}}) where {T,N} = IndexLinear()

Base.getindex(a::AbstractPooledDenseArray{T,N}, i::Int) where {T,N} = a.array[i]

Base.setindex!(a::AbstractPooledDenseArray{T,N}, v, i::Int) where {T,N} = a.array[i] = v


text:
# generate simple for loops for all dimensions

function generate_loops(Nmax, npts)
# Nmax is the maximum number of dimensions
# npts is the number of points in the stencil (used to figure out the name of
# the kernel to call)
# prefix is prepended to the file name

  str = ""
  for i=1:Nmax
    str_tmp = generateLoops(i, npts)
    str *= str_tmp
  end

  return str
end

function generateLoops(N, npts)
# generate a function that does N nested for loops

  np1 = N + 1
  str = ""
  str *= "function perf_hdindexing$npts{T}(u_i::AbstractArray{T, $np1}, u_ip1::AbstractArray{T, $np1})\n"

  indent = "  "

  # generate loop bounds
  for i=1:N
    varname = string("d", i, "min")
    varname2 = string("d", i, "max")
    idx = N - i + 1
    str *= indent*varname*" = 3\n"
    str *= indent*varname2*" = size(u_i, $idx) - 2\n"
  end

  str *= "\n"
  # generate loops
  for i=1:N
    varname = string("d", i)
    varname_min = string("d", i, "min")
    varname_max = string("d", i, "max")

    str *= indent*"for "*varname*" = "*varname_min*":"*varname_max*"\n"
    indent *= "  "
  end

  # generate body
  str_inner = ""
  str_inner *= indent*"idx = ("
  for i=1:N
    idx = N - i + 1  # indices are reversed
    varname = "d$idx"
    str_inner *= varname*", "
  end

  # remove trailing space and comma
  if N > 1
    str_inner = str_inner[1:end-2]
  end
  str_inner *= ")\n"
  str *= str_inner

  # call kernel
  str *= indent*"kernel$npts(idx, u_i, u_ip1)\n"

  # end loops
  for i=1:N
    indent = indent[1:end-2]
    str *= indent*"end\n"
  end

  str *= "\n"
  str *= indent*"return nothing\n"
  indent = indent[1:end-2]
  str *= "end\n"

  return str
end

#generate_loops(6, 5)


text:
"""
    SecretKey

Stores a secret key generated by a `KeyGenerator` instance.

See also: [`KeyGenerator`](@ref)
"""
mutable struct SecretKey <: SEALObject
  handle::Ptr{Cvoid}

  function SecretKey()
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:SecretKey_Create1, libsealc), Clong,
                   (Ref{Ptr{Cvoid}},),
                   handleref)
    @check_return_value retval
    return SecretKey(handleref[])
  end

  function SecretKey(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::SecretKey)
  if isallocated(object)
    @check_return_value ccall((:SecretKey_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function parms_id(key::SecretKey)
  parms_id = zeros(UInt64, 4)
  retval = ccall((:SecretKey_ParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 key, parms_id)
  @check_return_value retval
  return parms_id
end

function save!(buffer::DenseVector{UInt8}, length::Integer,
               compr_mode::ComprModeType.ComprModeTypeEnum, key::SecretKey)
  out_bytes = Ref{Int64}(0)
  retval = ccall((:SecretKey_Save, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}, UInt64, UInt8, Ref{Int64}),
                 key, buffer, length, compr_mode, out_bytes)
  @check_return_value retval
  return Int(out_bytes[])
end
function save!(buffer::DenseVector{UInt8}, length::Integer, key::SecretKey)
  return save!(buffer, length, ComprModeType.default, key)
end
function save!(buffer::DenseVector{UInt8}, key::SecretKey)
  return save!(buffer, length(buffer), key)
end

function save_size(compr_mode, key::SecretKey)
  result = Ref{Int64}(0)
  retval = ccall((:SecretKey_SaveSize, libsealc), Clong,
                 (Ptr{Cvoid}, UInt8, Ref{Int64}),
                 key, compr_mode, result)
  @check_return_value retval
  return Int(result[])
end
save_size(key::SecretKey) = save_size(ComprModeType.default, key)

function load!(key::SecretKey, context::SEALContext, buffer::DenseVector{UInt8}, length)
  in_bytes = Ref{Int64}(0)
  retval = ccall((:SecretKey_Load, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt8}, UInt64, Ref{Int64}),
                 key, context, buffer, length, in_bytes)
  @check_return_value retval
  return Int(in_bytes[])
end
function load!(key::SecretKey, context::SEALContext, buffer::DenseVector{UInt8})
  return load!(key, context, buffer, length(buffer))
end


text:
## Boundary conditions
# TODO: add class
# TODO: add constructors and other methods

## FunctionSpace class
@fenicsclass FiniteElement
export FiniteElement

## FunctionSpace class
@fenicsclass FunctionSpace
FunctionSpace(mesh::Mesh, family::String, degree::Int) = FunctionSpace(dolfin.FunctionSpace(mesh.pyobject, family, degree))
export FunctionSpace
# TODO: import operator(s), at least * (shall that be mapped to ⊗?)
const FAMILY_ARGYRIS = "ARG"
const FAMILY_ARNOLD_WINTHER = "AW"
const FAMILY_BREZZI_DOUGLAS_FORTIN_MARINI = "BDFM"
const FAMILY_BREZZI_DOUGLAS_MARINI = "BDM"
const FAMILY_BUBBLE = "B"
const FAMILY_CROUZEIX_RAVIART = "CR"
const FAMILY_DISCONTINUOUS_LAGRANGE = "DG"
const FAMILY_DISCONTINUOUS_RAVIART_THOMAS = "DRT"
const FAMILY_HERMITE = "HER"
const FAMILY_LAGRANGE = "CG"
const FAMILY_MARDAL_TAI_WINTHER = "MTW"
const FAMILY_MORLEY = "MOR"
const FAMILY_NEDELEC_1ST_KIND = "N1curl"
const FAMILY_NEDELEC_2ND_KIND = "N2curl"
const FAMILY_QUADRATURE = "Quadrature"
const FAMILY_RAVIART_THOMAS = "RT"
const FAMILY_REAL = "R"

mesh(V::FunctionSpace) = Mesh(fenicscall(V, :mesh))
ufl_element(V::FunctionSpace) = FiniteElement(fenicscall(V, :ufl_element))
dim(V::FunctionSpace) = fenicscall(V, :dim)::Int

export mesh, ufl_element, family, dim

## Expression class
@fenicsclass Expression
Expression(expr::String) = Expression(dolfin.Expression(expr))
Expression(expr::String, degree::Int) = Expression(dolfin.Expression(expr, degree=degree))
Expression(expr::String, element::FiniteElement) = Expression(dolfin.Expression(expr, element=element.pyobject))
Expression(expr::String, V::FunctionSpace) = Expression(expr, ufl_element(V))
export Expression
# TODO: add degree or element argument (make one necessary) (either degree=2, or V.ufl_element() or V directly)
# TODO: add optional domain arguments

## Expr class hierarchy (export only classes that are really needed in user code)
@fenicsclass Expr
@fenicsclass Argument Expr
@fenicsclass Constant Expr
TrialFunction(V::FunctionSpace) = Argument(dolfin.TrialFunction(V.pyobject))
TestFunction(V::FunctionSpace) = Argument(dolfin.TestFunction(V.pyobject))
Constant(x::Real) = Constant(dolfin.Constant(x, name="Constant($x)"))

export
  TrialFunction,
  TestFunction,
  Constant

## functions that operate on and return Expr's
nabla_grad(u::Expr) = Expr(dolfin.nabla_grad(u.pyobject))
inner(u::Expr, v::Expr) = Expr(dolfin.inner(u.pyobject, v.pyobject))
# TODO: add more functions (e.g. nabla_div, ...)
export
  nabla_grad, inner

## Measure class (export only concrete instances)
@fenicsclass Measure
dx = Measure(dolfin.dx)
ds = Measure(dolfin.ds)
export dx, ds

## Form class and operators that generate forms
@fenicsclass Form
*(expr::Expr, measure::Measure) = Form(measure.pyobject[:__rmul__](expr.pyobject) )
*(expr::Expr, expr2::Expr) = Expr(expr.pyobject[:__mul__](expr2.pyobject) )
# TODO: import more operators


text:
module EnKF

include("state.jl")

include("update.jl")

include("initial.jl")

include("inflation.jl")

include("system.jl")

include("stochasticEnKF.jl")

include("ETKF.jl")

end # module

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF


text:
import Base.Iterators: flatten

function _parseline(s::String)
    x = split(s, " <-> ")
    node = parse(Int, first(x))
    connected = parse.(Int, String.(split(x[2], ", ")))
    return Dict(node => connected)
end

function puzzle12(path::String=joinpath(@__DIR__, "..", "data/12.txt"))
    input = readlines(path)
    graph = reduce(merge, _parseline.(input))

    # start with 0
    nodes = [0]
    connected = nodes
    groups = Dict()

    while length(graph) > 0
        while true
            g = []
            for node in nodes
                append!(g, [get(graph, node, nothing)])
            end
            g = reduce(vcat, g)
            nodes = setdiff(g, connected)
            # nodes = setdiff(collect(flatten(get.(graph, nodes, nothing))), connected)
            if length(nodes) == 0
                break
            else
                for n in nodes
                    append!(connected, n)
                end
            end
        end
        groups[first(connected)] = connected
        for c in connected
            delete!(graph, c)
        end
        length(graph) == 0 ? break : nothing
        nodes = [first(keys(graph))]
        connected = nodes
    end

    partone = length(groups[0])
    parttwo = length(groups)

    return [partone, parttwo]
end


text:
#! /usr/bin/julia
#
#	text_read.jl
#
#					Jun/05/2018
# --------------------------------------------------------------------
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_in = ARGS[1]
println(fname_in)
#
dict_aa = text_to_dict_proc(fname_in)
dict_display_proc(dict_aa)
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
### Interactions are tested with the forest fire model
mutable struct Tree{T<:Integer} <: AbstractAgent
  id::T
  pos::Tuple{T, T}
  status::Bool  # true is green and false is burning
end

# we can put the model initiation in a function
function forest_initiation(;f, d, p, griddims, seed)
  Random.seed!(seed)

  space = Space(griddims, moore = true)

  properties = Dict(:f => f, :d => d, :p => p)
  forest = ABM(Tree, space; properties=properties, scheduler=random_activation)

  # create and add trees to each node with probability d, which determines the density of the forest
  for node in 1:nv(forest)
    pp = rand()
    if pp <= forest.properties[:d]
      tree = Tree(node, (1,1), true)
      add_agent!(tree, node, forest)
    end
  end
  return forest
end

function forest_step!(forest)
  for node in nodes(forest; by = :random)
    if length(forest.space.agent_positions[node]) == 0  # the cell is empty, maybe a tree grows here?
      p = rand()
      if p <= forest.properties[:p]
        bigest_id =  maximum(keys(forest.agents))
        treeid = bigest_id +1
        tree = Tree(treeid, (1,1), true)
        add_agent!(tree, node, forest)
      end
    else
      treeid = forest.space.agent_positions[node][1]  # id of the tree on this cell
      tree = id2agent(treeid, forest)  # the tree on this cell
      if tree.status == false  # if it is has been burning, remove it.
        kill_agent!(tree, forest)
      else
        f = rand()
        if f <= forest.properties[:f]  # the tree ignites on fire
          tree.status = false
        else  # if any neighbor is on fire, set this tree on fire too
          neighbor_cells = node_neighbors(tree, forest)
          for cell in neighbor_cells
            treeid = get_node_contents(cell, forest)
            if length(treeid) != 0  # the cell is not empty
              treen = id2agent(treeid[1], forest)
              if treen.status == false
                tree.status = false
                break
              end
            end
          end
        end
      end
    end
  end
end

@testset "Agent-Space interactions" begin

  model = forest_initiation(f=0.1, d=0.8, p=0.1, griddims=(20, 20), seed=2)  # forest fire model

  agent = model.agents[1]
  move_agent!(agent, (3,4), model)  # node number 63
  @test agent.pos == (3,4)
  @test agent.id in model.space.agent_positions[63]

  agent = model.agents[2]
  move_agent!(agent, 83, model)  # pos (3,5)
  @test agent.pos == (3,5)
  @test agent.id in model.space.agent_positions[83]

  new_pos = move_agent!(agent, model)
  @test agent.id in model.space.agent_positions[coord2vertex(new_pos, model)]

  add_agent!(agent, (2,9), model)
  @test agent.pos == (2,9)
  @test agent.id in model.space.agent_positions[coord2vertex((2,9), model)]
  @test agent.id in model.space.agent_positions[coord2vertex(new_pos, model)]

  model1 = ABM(Agent1, Space((3,3)))
  add_agent!(1, model1)
  @test model1.agents[1].pos == (1, 1)
  add_agent!((2,1), model1)
  @test model1.agents[2].pos == (2, 1)

  model2 = ABM(Agent4, Space((3,3)))
  add_agent!(1, model2, 3)
  @test model2.agents[1].pos == (1,1)
  @test 1 in model2.space.agent_positions[1]
  add_agent!((2,1), model2, 2)
  @test model2.agents[2].pos == (2,1)
  @test 2 in model2.space.agent_positions[2]
  ag = add_agent!(model2, 12)
  @test ag.id in model2.space.agent_positions[coord2vertex(ag.pos, model2)]

  @test agent.id in get_node_contents(agent, model)

  ii = model.agents[length(model.agents)]
  @test id2agent(ii.id, model) == model.agents[ii.id]

  agent = model.agents[1]
  agent_pos = coord2vertex(agent.pos, model)
  kill_agent!(agent, model)
  @test_throws KeyError id2agent(1, model)
  @test !in(1, Agents.agent_positions(model)[agent_pos])
end


text:
@testset "275.h-index.jl" begin
    @test h_index_ii([0,1,3,5,6]) == 3
end


text:
function show_dag_path(d::DAG, path::Vector{Symbol})
  local str
  for (ind, sym) in enumerate(path)
    if ind == 1
      str = ":$sym"
    else
      cond = d.a[path[ind-1], path[ind]] == 1
      str *= cond ? " \u21d0 :$(path[ind])" : " \u21d2 :$(path[ind])"
    end
  end
  str
end

export
  show_dag_path


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the beta angle of a satellite.
#
#   The algorithm was based on:
#
#       Mortari, D., Wilkins, M. P., and Bruccoleri, C.  On Sun-Synchronous
#       Orbits and Associated Constellations
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export beta_angle

"""
    beta_angle(JD₀::Number, a::Number, e::Number, i::Number, RAAN::Number, Δt::Integer, pert::Symbol = :J2)

Compute the beta angle of an orbit with semi-major axis `a` [m], eccentricity
`e`, inclination `i` [rad], and initial right ascension of the ascending node
`RAAN` [rad]. The orbit epoch, which is also the day in which the analysis will
begin, is `JD₀` [Julian Day]. The analysis will be performed for each day during
`Δt` days.

The argument `pert` can be used to select the perturbation terms that must be
used when propagating the right ascencion of the ascending node. The possible
values are:

* `:J0`: Consider a Keplerian orbit.
* `:J2`: Consider the perturbation terms up to J2.
* `:J4`: Consider the perturbation terms J2, J4, and J2².

If `pert` is omitted, then it defaults to `:J2`.

# Returns

An array with two columns. The first one contains the days of the analysis and
the second one contains the beta angle [rad] for each day.

"""
function beta_angle(JD₀::Number, a::Number, e::Number, i::Number, RAAN::Number,
                    Δt::Integer, pert::Symbol = :J2)

    # Vector of the days in which the beta angle will be computed.
    days = collect(0:1:Δt-1)

    # Output vector.
    β = Vector{Float64}(undef,Δt)

    # RAAN rotation rate [rad/day].
    δΩ = 86400dRAAN(a, e, i, pert)

    # Loop
    @inbounds @views for d in days
        # Compute the RAAN at the day d.
        Ω = RAAN + δΩ*d

        # Compute the versor N represented in the Inertial ref. frame.
        Dio = angle_to_dcm(-i, -Ω, 0, :XZX)
        N_i = Dio[:,3]  # -> Dio*[0;0;1]

        # Compute the Sun position at noon (UT) represented in the Inertial ref.
        # frame.
        S_i = sun_position_i(JD₀+d)
        S_i = S_i/norm(S_i)

        # Compute the beta angle, which is the angle between the Sun vector and
        # the orbit plane.
        β[d+1] = 90 - acosd(dot(N_i,S_i))
    end

    return [days β]
end


text:
struct GeoPlot end

# We only look at arrays with Lat/Lon here.
# Otherwise they fall back to DimensionalData.jl recipes
@recipe function f(A::AbstractGeoArray)
    A = GeoArray(A)
    if all(hasdim(A, (Lat(), Lon())))
        # Heatmap or multiple heatmaps. Use GD recipes.
        GeoPlot(), prepare(A)
    else
        # This is not a Lat/Lon heatmap. Fall back to DD recipes after reprojecting
        da = A |> GeoArray |> a -> DimArray(a; dims=_maybe_mapped(dims(a)))
        DimensionalData.DimensionalPlot(), da
    end
end

# Plot 3d arrays as multiple tiled plots
@recipe function f(::GeoPlot, A::GeoArray{T,3,<:Tuple{<:Lat,<:Lon,D}}) where {T,D}
    nplots = size(A, 3)
    if nplots > 1
        :layout --> nplots
        :xguide --> permutedims(string.(val(dims(A, D))))
        for i in 1:nplots
            @series begin
                seriestype := :heatmap
                aspect_ratio := 1
                subplot := i
                slice = A[:, :, i]
                lats, lons = map(prepare, dims(slice))
                lons, lats, parent(slice)
            end
        end
    else
        GeoPlot(), A[:, :, 1]
    end
end

# # Permute for correct Lat/Lon order
@recipe function f(::GeoPlot, A::GeoArray{T,3}) where {T}
    GeoPlot(), permutedims(A, (Lat, Lon, Dimension))
end

# # Plot a sinlge 2d map
@recipe function f(::GeoPlot, A::GeoArray{T,2,<:Tuple{<:Lat,<:Lon}}) where T
    :seriestype --> :heatmap
    :aspect_ratio --> 1
    :colorbar_title --> name(A)
    :title --> DD._refdims_title(A)
    lats, lons = map(prepare, dims(A))
    lons, lats, parent(A)
end

# # Permute for correct Lat/Lon order
@recipe function f(::GeoPlot, A::GeoArray{T,2,<:Tuple{<:Lon,<:Lat}}) where T
    GeoPlot(), permutedims(A)
end

# Plots heatmaps pixels are centered.
# So we should center, and use the projected value.
prepare(d::Dimension) = shiftindexloci(Center(), d) |> _maybe_mapped |> index

# Convert arrays to a consistent missing value and Forward array order
prepare(A::AbstractGeoArray) =
    _maybe_replace_missing(A) |>
    A -> reorder(A, ForwardIndex) |>
    A -> reorder(A, ForwardRelation)

_maybe_replace_missing(A::AbstractArray{<:AbstractFloat}) = replace_missing(A, eltype(A)(NaN))
_maybe_replace_missing(A) = A

_maybe_mapped(dims::Tuple) = map(_maybe_mapped, dims)
_maybe_mapped(dim::Dimension) = _maybe_mapped(mode(dim), dim)
_maybe_mapped(mode::IndexMode, dim::Dimension) = dim
_maybe_mapped(mode::Projected, dim::Dimension) = _maybe_mapped(mappedcrs(mode), dim)
_maybe_mapped(::Nothing, dim::Dimension) = dim
_maybe_mapped(::GeoFormat, dim::Dimension) = convertmode(Mapped, dim)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestKVStore
using MXNet
using Test

using ..Main: rand_dims

SHAPE = rand_dims()
KEYS  = [5,7,11]

function init_kv()
  kv = mx.KVStore()
  mx.init!(kv, 3, mx.zeros(SHAPE))

  vals = [mx.zeros(SHAPE) for k in KEYS]
  mx.init!(kv, KEYS, vals)
  kv
end

function test_kv_basic()
  @info("KVStore::basic")

  kv = init_kv()
  @test mx.get_type(kv) == :local
  @test mx.get_rank(kv) == 0
  @test mx.get_num_workers(kv) == 1
end

function test_single_kv_pair()
  @info("KVStore::single")

  kv = init_kv()
  mx.push!(kv, 3, mx.ones(SHAPE))
  val = NDArray(undef, SHAPE)
  mx.pull!(kv, 3, val)
  @test maximum(abs.(copy(val) .- 1)) == 0
end

function test_aggregator()
  @info("KVStore::aggregator")

  kv = init_kv()

  num_devs = 4
  devs = [mx.Context(mx.CPU, i) for i=0:num_devs-1]
  vals = [mx.ones(SHAPE, dev) for dev in devs]

  mx.push!(kv, 3, vals)
  mx.pull!(kv, 3, vals)
  for v in vals
    @test maximum(abs.(copy(v)) .- num_devs) == 0
  end

  # list
  vals = [mx.NDArray[mx.ones(SHAPE, dev)*2 for dev in devs] for k in KEYS]
  mx.push!(kv, KEYS, vals)
  mx.pull!(kv, KEYS, vals)

  for vv in vals
    for v in vv
      @test maximum(abs.(copy(v)) .- 2 * num_devs) == 0
    end
  end
end

function check_setupdater!(f)
  kv = KVStore(:local)
  setupdater!(kv, f)

  A = Float32[1, 2, 3, 4]
  B = Float32[.5, .6, .7, .8]
  x = NDArray(A)
  Δ = NDArray(B)
  init!(kv, 42, x)
  push!(kv, 42, Δ)
  pull!(kv, 42, x)

  @test copy(x) ≈ A + 2B
end  # function check_setupdater!

function test_setupdater!()
  @info("KVStore::setupdater!")

  f(key, Δ, x) = @mx.inplace x += 2Δ
  g(key, Δ, x) = (x[:] += 2Δ)

  check_setupdater!(f)
  check_setupdater!(g)
end  # test_setupdater!

@testset "KVStore Test" begin
  test_kv_basic()
  test_single_kv_pair()
  test_aggregator()
  test_setupdater!()
end

end


text:
"""
Overview: parametrisation structure management
"""

"""
> mutable struct Parametrisation
- nc : number of monomial combinations
- W : mappings [Υ, Ψ]
- f : reduced dynamics
- nlr : nonlinear right hand side
- comb : combinations the system is solved for
- conj : conjugacy relations
- cmap : mapping index combinations
"""
mutable struct Parametrisation
  nc::Int64
  W::Matrix{ComplexF64}
  f::Matrix{ComplexF64}
  nlr::Matrix{ComplexF64}
  comb::Matrix{Int64}
  conj::Array{Int64}
  cmap::Array{Int64}
  #
  Parametrisation() = new()
  #
end

"""
> initialize_parametrisation!(param,order,ndofs,neq)
It initializes a generic order of the asymptotic expansion
- param : parametrisation
- order : order
- ndofs : number of degrees of freedom of the reduced model
- neq   : number of equations
"""
function initialize_parametrisation!(param::Parametrisation,p::Int64,
                                     ndofs::Int64,neq::Int64,veps=0)
  #
  nperm = ndofs^p
  nc  = floor(Int64,compute_ncomb(p,ndofs))
  param.nc = nc
  if (p>0)
    comb = zeros(Int64,p,nc)
  else
    comb = zeros(Int64,1,nc)
  end
  conj = zeros(Int64,nc)
  cmap = zeros(Int64,nperm)
  #
  if (p>1)
    #
    ith_c = 1
    combs = zeros(Int64,p)
    sorted_combs = zeros(Int64,p)
    #
    ith_c = fill_comb!(comb,cmap,ndofs,p,combs,ith_c,1)
    generate_pmap!(comb,cmap,ndofs,p,combs,sorted_combs,1)
    generate_conj!(conj,comb,cmap,ndofs,p,combs,nc)
    #
  elseif (p==1)
    nm = Int(ndofs/2)
    for i = 1:nm
      conj[i] = i
      conj[i+nm] = -i
    end
    for i = 1:ndofs
      comb[p,i] = i
      cmap[i] = i
    end
  else
    conj[1] = 1
    comb[1,1] = 0
    cmap[1] = 1
  end
  #
  if (veps==1 && p>0)
    [conj[i]=i for i in 1:p]
  end
  #
  param.comb = comb
  param.conj = conj
  param.cmap = cmap
  #
  param.W   = zeros(ComplexF64,neq*2,nc)
  param.nlr = zeros(ComplexF64,neq,nc)
  param.f   = zeros(ComplexF64,ndofs,nc)
  #
  return nothing
end

"""
> generate_conj!(conj,comb,cmap,ndofs,p,combs,nc)
It maps entries to their conjugate
- conj : conjugate
- comb : matrix that stores the powers of the monomials
- cmap : total combinations map
- ndofs : number of degrees of freedom
- p : order of the asymptotic expansion
- combs : combinations
- nc : number of combinations
"""
function generate_conj!(conj,comb,cmap,ndofs,p,combs,nc)
  t = 1
  nm = Int(ndofs/2)
  for i = 1:nc
    #
    if (conj[i]!=0)
      continue
    end
    #
    for j = 1:p
      combs[j] = comb[j,i]
    end
    #
    conj[i] = t
    #
    for j = 1:p
      if (combs[j]>nm)
        combs[j] -= nm
      else
        combs[j] += nm
      end
    end
    sort!(combs)
    #
    pos1 = find_index_global(combs,ndofs,p)
    #
    cm = cmap[pos1]
    if (cm>0 && conj[cm]==0)
      conj[cm] = -i
    end
    t += 1
  end
  #
  return nothing
end

"""
> generate_pmap!(comb::Matrix{Int64},cmap::Array{Int64},ndofs::Int64,p::Int64,combs::Array{Int64},sorted_combs::Array{Int64},counter::Int64)
It computes combinations maps to account for indexes permutations.
- comb : matrix that stores the powers of the monomials
- cmap : total combinations map
- ndofs : number of degrees of freedom
- p : order of the asymptotic expansion
- combs : combinations
- sorted_combs : sorted combinations
- counter : recursion depth
"""
function generate_pmap!(comb::Matrix{Int64},cmap::Array{Int64},ndofs::Int64,p::Int64,
                        combs::Array{Int64},sorted_combs::Array{Int64},counter::Int64)
  for i = 1:ndofs
    combs[counter] = i
    if (counter<p)
      counter += 1
      generate_pmap!(comb,cmap,ndofs,p,combs,sorted_combs,counter)
      counter -= 1
    else
      pos1 = find_index_global(combs,ndofs,p)
      sorted_combs = sort(combs)
      pos2 = find_index_global(sorted_combs,ndofs,p)
      if !(pos1==pos2)
        cmap[pos1] = -cmap[pos2]
      end
    end
  end
end

"""
> fill_comb!(comb::Matrix{Int64},cmap::Array{Int64},ndofs::Int64,p::Int64,combs::Array{Int64},ith_c::Int64,counter::Int64)
It fills the comb matrix with the associated exponentials
- comb : matrix that stores the powers of the monomials
- cmap : total combinations map
- ndofs : number of degrees of freedom
- p : order of the asymptotic expansion
- combs : combinations
- ith_c : ith-combination
- counter : recursion depth
"""
function fill_comb!(comb::Matrix{Int64},cmap::Array{Int64},ndofs::Int64,
                    p::Int64,combs::Array{Int64},ith_c::Int64,counter::Int64)
  #
  if (counter==1)
    for i = 1:ndofs
      combs[counter] = i
      counter += 1
      ith_c = fill_comb!(comb,cmap,ndofs,p,combs,ith_c,counter)
      counter -= 1
    end
  else
    for i = combs[counter-1]:ndofs
      combs[counter] = i
      if (counter<p)
        counter += 1
        ith_c = fill_comb!(comb,cmap,ndofs,p,combs,ith_c,counter)
        counter -= 1
      else
        for j = 1:p
          comb[j,ith_c] = combs[j]
        end
        pos = find_index_global(combs,ndofs,p)
        cmap[pos] = ith_c
        ith_c += 1
      end
    end
  end
  return ith_c
end

"""
> compute_ncomb(param,order,ndofs,neq)
It computes the monomials combinations that need to be computed
- p : order of the expansion
- ndofs : number of degrees of freedom of the reduced model
"""
function compute_ncomb(p::Int64,ndofs::Int64)
  #
  tmp = p+ndofs-1
  nc = 1
  for i = 1:tmp
    nc = nc*i
    if (i<=p)
      nc /= i
    end
  end
  for i = 1:ndofs-1
    nc /= i
  end
  return nc
  #
end

function free_parametrisation!(param::Parametrisation)
  #
  param.W = zeros(ComplexF64,(1,1))
  param.f = zeros(ComplexF64,(1,1))
  param.nlr = zeros(ComplexF64,(1,1))
  param.comb = zeros(Int64,(1,1))
  param.conj = zeros(Int64,1)
  param.cmap = zeros(Int64,1)
  #
  GC.gc()
  #
  return nothing

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LinearPartition"
version = v"1.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/LinearFold/LinearPartition/archive/refs/tags/v$(version.major).$(version.minor).tar.gz",
                  "4fdea96f7ffbd4804d9308ddb46db5f96d1abc4b7bd737725f9bedcae3c88178"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/LinearPartition-*/

make -j${nproc} CC=${CXX}
mkdir -p ${bindir}
for prg in linearpartition_{c,v}; do
    install -Dvm 755 "bin/${prg}" "${bindir}/${prg}${exeext}"
done

install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("linearpartition_c", :linearpartition_c),
    ExecutableProduct("linearpartition_v", :linearpartition_v),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")


text:
module BlockDiagonalFactors

using SparseArrays, SuiteSparse, LinearAlgebra

struct SparseBlockFactors{Tv}
    factors::Vector
    indices::Vector{<:Int}
    m::Int
    n::Int
end
struct BlockFactors{Tv}
    factors::Vector
    indices::Vector{<:Int}
    m::Int
    n::Int
end
export BlockFactors, SparseBlockFactors

# Factorization functions
for f in (:lu, :qr, :cholesky, :factorize), T in (:Float64, :(Complex{Float64}))
    @eval begin
        import LinearAlgebra: $f
        """
        $($f)(As::Vector{SparseMatrixCSC{$($T),Int}}, I::Vector{Int})

        Creates a block-diagonal (lazy) array of factors.
        Invokes `$($f)` on each matrix in the array of matrices `As` and stores them along with the indices `I`.
        """
        function $f(As::Vector{SparseMatrixCSC{$T,Int}}, I::Vector{Int})
            m = sum(i -> As[i].m, I)
            n = sum(i -> As[i].n, I)
            return SparseBlockFactors{$T}(map($f, As), I, m, n)
        end
        """
        $($f)(As::Vector{Array{$($T),2}}, I::Vector{Int})

        Creates a block-diagonal (lazy) array of factors.
        Invokes `$($f)` on each matrix in the array of matrices `As` and stores them along with the indices `I`.
        """
        function $f(As::Vector{Array{$T,2}}, I::Vector{Int})
            m = sum(i -> size(As[i],1), I)
            n = sum(i -> size(As[i],2), I)
            return BlockFactors{$T}(map($f, As), I, m, n)
        end
        export $f
    end
end

# TODO replace with lazy version? not sure how to do this
for f in (:adjoint, :transpose), TF in (:SparseBlockFactors, :BlockFactors), T in (:Float64, :(Complex{Float64}))
    @eval begin
        import Base: $f
        """
        $($f)(BDF::$($T))

        Invokes `$($f)` on all the factors in `BDF` and returns them into a new `$($T)` object.
        """
        $f(BDF::$TF{$T}) = $TF{$T}(map($f, BDF.factors), BDF.indices, BDF.m, BDF.n)
        export $f
    end
end

import Base.\
"""
    \\(BDF::SparseBlockFactors, y::AbstractVecOrMat)

Backsubstitution for `SparseBlockFactors`.
"""
function \(BDF::SparseBlockFactors{T}, y::AbstractVector{S}) where {T,S}
    x = Array{promote_type(T,S)}(undef, BDF.n)
    x_idx, y_idx = 0:0, 0:0
    my = size(y)[1]
    for i in BDF.indices
        x_idx = x_idx.stop .+ (1:size(BDF.factors[i],2))
        y_idx = y_idx.stop .+ (1:size(BDF.factors[i],1))
        x[x_idx] .= BDF.factors[i] \ y[y_idx]
    end
    return x
end
function \(BDF::SparseBlockFactors{T}, y::AbstractVector{T}) where T
    x = Array{T}(undef, BDF.n)
    x_idx, y_idx = 0:0, 0:0
    my = size(y)[1]
    for i in BDF.indices
        x_idx = x_idx.stop .+ (1:size(BDF.factors[i],2))
        y_idx = y_idx.stop .+ (1:size(BDF.factors[i],1))
        x[x_idx] .= BDF.factors[i] \ y[y_idx]
    end
    return x
end
"""
    \\(BDF::BlockFactors, y::AbstractVecOrMat)

Backsubstitution for `BlockFactors`.
"""
function \(BDF::BlockFactors{T}, y::AbstractVector{S}) where {T,S}
    x = Array{promote_type(T,S)}(undef, BDF.n)
    x_idx, y_idx = 0:0, 0:0
    my = size(y)[1]
    for i in BDF.indices
        x_idx = x_idx.stop .+ (1:size(BDF.factors[i],2))
        y_idx = y_idx.stop .+ (1:size(BDF.factors[i],1))
        x[x_idx] .= BDF.factors[i] \ y[y_idx]
    end
    return x
end
function \(BDF::BlockFactors{T}, y::AbstractVector{T}) where T
    x = Array{T}(undef, BDF.n)
    x_idx, y_idx = 0:0, 0:0
    my = size(y)[1]
    for i in BDF.indices
        x_idx = x_idx.stop .+ (1:size(BDF.factors[i],2))
        y_idx = y_idx.stop .+ (1:size(BDF.factors[i],1))
        x[x_idx] .= BDF.factors[i] \ y[y_idx]
    end
    return x
end
export \

import LinearAlgebra.ldiv!
"""
    ldiv!(BDF::SparseBlockFactors, y::AbstractVecOrMat)

Backsubstitution for `SparseBlockFactors`. For square blocks only.
"""
function ldiv!(BDF::SparseBlockFactors{T}, y::AbstractVector{T}) where T
    idx = 0:0
    for i in BDF.indices
        idx = idx.stop .+ (1:BDF.factors[i].n)
        y[idx] .= ldiv!(BDF.factors[i], y[idx])
    end
    return y
end
"""
    ldiv!(BDF::BlockFactors, y::AbstractVecOrMat)

Backsubstitution for `BlockFactors`. For square blocks only.
"""
function ldiv!(BDF::BlockFactors{T}, y::AbstractVector{T}) where T
    idx = 0:0
    for i in BDF.indices
        idx = idx.stop .+ (1:size(BDF.factors[i],2))
        y[idx] .= ldiv!(BDF.factors[i], y[idx])
    end
    return y
end
"""
    ldiv!(x::AbstractVecOrMat, BDF::SparseBlockFactors, y::AbstractVecOrMat)

Backsubstitution for `SparseBlockFactors`. For square blocks only.
"""
function ldiv!(x::AbstractVecOrMat{T}, BDF::SparseBlockFactors{T}, y::AbstractVector{T}) where T
    idx = 0:0
    for i in BDF.indices
        idx = idx.stop .+ (1:BDF.factors[i].n)
        x[idx] .= ldiv!(x[idx], BDF.factors[i], y[idx])
    end
    return x
end
"""
    ldiv!(x::AbstractVecOrMat, BDF::BlockFactors, y::AbstractVecOrMat)

Backsubstitution for `BlockFactors`. For square blocks only.
"""
function ldiv!(x::AbstractVecOrMat{T}, BDF::BlockFactors{T}, y::AbstractVector{T}) where T
    idx = 0:0
    for i in BDF.indices
        idx = idx.stop .+ (1:size(BDF.factors[i],2))
        x[idx] .= ldiv!(x[idx], BDF.factors[i], y[idx])
    end
    return x
end
export ldiv!

end # module


text:
using QuantumCumulants
using SymbolicUtils
using Test

@testset "meanfield" begin

hf = FockSpace(:cavity)
ha = NLevelSpace(:atom,(:g,:e))
h = hf⊗ha

a = Destroy(h,:a)
σ = Transition(h,:σ,:g,:e)
σee = Transition(h,:σ,:e,:e,2)

# JC
(ωc,ωa,g) = (1.1341,0.4321,2.15013)
H = ωc*a'*a + ωa*σ'*σ + g*(a'*σ + σ'*a)

da = commutator(im*H,a)
@test iszero(simplify(da - (-im*ωc*a + (-im*g)*σ)))
ds = commutator(im*H,σ)
@test iszero(simplify(ds - ((-im*g)*a + (-im*ωa)*σ + (2im*g)*a*σee)))
dn = commutator(im*H,a'*a)
@test iszero(simplify(dn - ((-im*g)*a'*σ + (im*g)*a*σ')))

he = meanfield([a,σ,a'*a],H)
@test iszero(simplify(he.operator_equations[1].rhs - (da)))
@test iszero(simplify(he.operator_equations[2].rhs - (ds)))
@test iszero(simplify(he.operator_equations[3].rhs - (dn)))

# Lossy JC
κ,γ = 3.333,0.1313131313
J = [a,σ]
he_diss = meanfield([a,σ,σ'*σ],H,J;rates=[κ,γ])

@test iszero(simplify(he_diss.operator_equations[1].rhs - ((-im*ωc - 0.5κ)*a + (-im*g)*σ)))
@test iszero(simplify(he_diss.operator_equations[2].rhs - ((-im*g)*a + (-im*ωa - 0.5γ)*σ + (2im*g)*a*σee)))
@test iszero(simplify(he_diss.operator_equations[3].rhs - ((-γ)*σee + (im*g)*a'*σ + (-im*g)*a*σ')))

# Single-atom laser
ν = 3.44444444
J = [a,σ,σ']
he_laser = meanfield([a'*a,σ'*σ,a'*σ],H,J;rates=[κ,γ,ν])

@test iszero(simplify(he_laser.operator_equations[1].rhs - ((-κ)*a'*a + (-im*g)*a'*σ + (im*g)*a*σ')))
@test iszero(simplify(he_laser.operator_equations[2].rhs - (ν + (-ν - γ)*σee + (im*g)*a'*σ + (-im*g)*a*σ')))
ex = (im*g)*σee + (-im*g)*a'*a + (im*(ωc - ωa) - 0.5*(κ + γ + ν))*a'*σ + (2im*g)*a'*a*σee
@test iszero(simplify(he_laser.operator_equations[3].rhs - ex))

# collective decay
N=2
@cnumbers G11 G12 G21 G22 δ
h = ⊗([NLevelSpace(Symbol(:atom,i),2) for i=1:N]...)
σ_(i,j,k) = Transition(h,Symbol("σ__{$k}"),i,j,k)
H = δ*(σ_(2,2,1)+σ_(2,2,2))
J = [[σ_(1,2,1), σ_(1,2,2)]]
rates = [[G11 G12; G21 G22]]
ops_0 = [σ_(1,2,1)]
eqs1 = meanfield(ops_0,H,J; rates=rates)
eqs_c1 = complete(eqs1)
eqs2 = meanfield(ops_0,H,J[1]; rates=rates[1])
eqs_c2 = complete(eqs2)
@test length(eqs_c2) == 4

JumpOp = Transition{ProductSpace{Vector{NLevelSpace{Symbol, UnitRange{Int64}, Int64}}}, Symbol, Int64, Int64}[]
JumpOpConj = Transition{ProductSpace{Vector{NLevelSpace{Symbol, UnitRange{Int64}, Int64}}}, Symbol, Int64, Int64}[]
for i=1:N
    for j=1:N
        push!(JumpOp,σ_(1,2,i))
        push!(JumpOpConj,σ_(2,1,j))
    end
end
rates=[G11, G21, G12, G22]
eqs3 = meanfield(ops_0,H,JumpOp;Jdagger=JumpOpConj, rates=rates)
eqs_c3 = complete(eqs3)

ops4_0 = QuantumCumulants.undo_average.(eqs_c1.states)
eqs4 = meanfield(ops4_0,H,JumpOp;Jdagger=JumpOpConj, rates=rates)
eqs_c4 = complete(eqs4)
@test eqs_c1.equations == eqs_c2.equations == eqs_c3.equations == eqs_c4.equations # jumps_dagger are different

end # testset


text:
module Keyboards

using Parameters

using Unitful
import Unitful: Length

using Circuits
using Circuits.TikZ
using TikzPictures
const preamble = read(joinpath(dirname(@__FILE__), "preamble.tex"), String)

using PCBs

include("design.jl")
include("circuits.jl")

end # module


text:
using FirebirdDB
using Base.Test

include("firebirddb_tests.jl")


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 1f43297c-5e20-11eb-1ce3-5b75c0b5c368
begin
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
	using StatsBase
end

# ╔═╡ 5818518e-5e20-11eb-3a74-c3f0b78d8098
md"""
## Left Right Environment
"""

# ╔═╡ 6fe63566-5e20-11eb-056c-1f513ae5d161
begin
	Base.@kwdef mutable struct LeftRightEnv <: AbstractEnv
		reward::Float64 = 0.
		current_state::Int = 1
	end

	RLBase.state_space(env::LeftRightEnv) = Base.OneTo(2)
	RLBase.action_space(env::LeftRightEnv) = Base.OneTo(2)

	function (env::LeftRightEnv)(a::Int)
		if a == 2
			env.reward = 0.
			env.current_state = 2
		else
			s = sample(Weights([0.9, 0.1], 1.0))
			if s == 1
				env.reward = 0.
				env.current_state = 1
			else
				env.reward = 1.0
				env.current_state = 2
			end
		end
	end

	function RLBase.reset!(env::LeftRightEnv)
		env.current_state = 1
		env.reward = 0.
	end

	RLBase.reward(env::LeftRightEnv) = env.reward
	RLBase.is_terminated(env::LeftRightEnv) = env.current_state == 2
	RLBase.state(env::LeftRightEnv) = env.current_state
end

# ╔═╡ 54f876d4-5e21-11eb-14be-67265069c8ae
world = LeftRightEnv()

# ╔═╡ 62022db4-5e21-11eb-220d-d163a3aa5c81
begin
	ns = length(state_space(world))
	na = length(action_space(world))
end

# ╔═╡ 597bbcd2-5e21-11eb-10c9-7db2f46cfffb
 π_t = VBasedPolicy(
	learner=MonteCarloLearner(
		approximator=(
			TabularVApproximator(;n_state=ns, opt=Descent(1.0)), # V
			TabularVApproximator(;n_state=ns, opt=InvDecay(1.0)) # Returns
			),
		kind=FIRST_VISIT,
		sampling=ORDINARY_IMPORTANCE_SAMPLING,
		γ=1.0
		),
	mapping= (env, V) -> 1
)

# ╔═╡ 9cdef3ba-5e21-11eb-1cfd-fb1abbf1d608
# A little ad-hoc here
RLBase.prob(::typeof(π_t), s, a) = a == 1 ? 1.0 : 0.

# ╔═╡ 88b48522-5e20-11eb-00df-65834ec124b2
begin
	struct CollectValue <: AbstractHook
		values::Vector{Float64}
		CollectValue() = new([])
	end
	(f::CollectValue)(::PostEpisodeStage, agent, env) = push!(f.values, agent.policy.π_target.learner.approximator[2](1))
end

# ╔═╡ 4e0a26b0-5e21-11eb-01bf-85db236b9bf8
begin
	p = plot()
	for _ in 1:10
		fill!(π_t.learner.approximator[1].table, 0.)
		fill!(π_t.learner.approximator[2].table, 0.)
		empty!(π_t.learner.approximator[2].optimizer.state)
		agent = Agent(
			policy=OffPolicy(
				π_t,
				RandomPolicy(1:na)
				),
			trajectory=VectorWSARTTrajectory()
		)
		hook = CollectValue()
		run(agent, world, StopAfterEpisode(100_000, is_show_progress=true),hook)
		plot!(p, hook.values, xscale = :log10)
	end
	p
end

# ╔═╡ Cell order:
# ╠═1f43297c-5e20-11eb-1ce3-5b75c0b5c368
# ╟─5818518e-5e20-11eb-3a74-c3f0b78d8098
# ╠═6fe63566-5e20-11eb-056c-1f513ae5d161
# ╠═54f876d4-5e21-11eb-14be-67265069c8ae
# ╠═62022db4-5e21-11eb-220d-d163a3aa5c81
# ╠═597bbcd2-5e21-11eb-10c9-7db2f46cfffb
# ╠═9cdef3ba-5e21-11eb-1cfd-fb1abbf1d608
# ╠═88b48522-5e20-11eb-00df-65834ec124b2
# ╠═4e0a26b0-5e21-11eb-01bf-85db236b9bf8


text:
#TODO implement proper soil module
@system Soil begin
    T_soil => 10 ~ track(u"°C")
    WP_leaf: leaf_water_potential => 0 ~ track(u"MPa") # pressure - leaf water potential MPa...
    total_root_weight => 0 ~ track(u"g")
end

@system SoilStub begin
    soil ~ hold

    WP_leaf(soil.WP_leaf) ~ track(u"MPa")
end

# Base.show(io::IO, s::Soil) = print(io, "$(s.T_soil)\n$(s.WP_leaf)\n$(s.total_root_weight)")


text:
module RAIGenerated
using InteractiveUtils

include("utilities.jl")

"""
    @rAI_generated function f(x) :(x+1) end

Implements the same functionality as @generated, but using `@eval`, so that it is guaranteed
to generate a new function for given type parameters exactly once.
"""
macro rAI_generated(f)
    if isa(f, Expr) && (f.head === :function || Base.is_short_function_def(f))
        body = f.args[2]
        lno = body.args[1]
        return gen_eval_f(f)
    else
        error("invalid syntax; @generated must be used with a function definition")
    end
end

function raw_argnames(args::Array)
    tmpcount = 0
    out = []
    for a in args
        name = raw_argname(a)
        if name == nothing
            tmpcount += 1
            name = Symbol("_$tmpcount")
        end
        push!(out, name)
    end
    out
end
raw_argname(x::Symbol) = x
function raw_argname(e::Expr)
    @assert e.head == Symbol("::")
    return length(e.args) == 2 ? e.args[1] : nothing
end

function gen_eval_f(funcexpr::Expr)
    dump(funcexpr)
    name = funcexpr.args[1].args[1]
    escname = esc(name)
    args = funcexpr.args[1].args[2:end]
    types = argtypes(args)
    names = raw_argnames(args)
    esc_argnames = argnames(args)
    function_params = paramexprs(args)
    escargs = [esc(a) for a in args]
    body = esc(funcexpr.args[2])
    quote
        function funchelper($(escargs...))
            $body
        end
        function $escname($(function_params...))
            names = $names
            types = [$((:typeof(x) for x in params)...)]
            generated_body = funchelper(types...)
            typedargs = [Expr(Symbol("::"), names[i], types[i]) for i in 1:$(length(args))]
            #@eval $finalcall = $(Expr(:$, :val))
            @eval $name($(Expr(:$, :(typedargs...)))) = $(Expr(:$, :generated_body))
            return @eval $name($([Expr(:$, n) for n in names]...))
        end
    end
end
gen_eval_f(:(f(x::Int64) = 3))

# EXAMPLE:
# This is what the macro expands to:
function bar(x)
    names = [:x]
    types = [typeof(x)]
    gen_body = bar_helper(types...)
    # Real @generated creates method instances, not methods
    #  Consequence: allows you to redefine body of bar() in the future.
    # Also: slower! Dynamic dispatch is slow, and gets slower w/ more methods
    @eval bar(x::types...) = $(gen_body)
    return @eval bar(x)
end

@rAI_generated function bar(x)
    if x <: Integer
        return :(x ^ 2)
    else
        return :(x)
    end
end

# Generated functions in general:
#  - dynamic dispatch, more often, and with more methodinstances
#  - Forced to have one methodinstances per unique params type tuple
#  - Forced to know the runtime types always to have statically inferrable code

@show methods(bar)
@show bar(3.0)
@show methods(bar)
@show bar(3.0)
@show bar(3)
@code_typed bar(3.0)

# Compare with the built-in @generated, which doesn't add new methods:
# @generated function baz(x)
#     if x <: Integer
#         return :(x ^ 2)
#     else
#         return :(x)
#     end
# end
#
# @show methods(baz)
# baz(4)

end


text:
function Source{T1<:AbstractFloat,T2<:Complex}(isz::Array{Int,1},isx::Array{Int,1},ot::Array{T1,1},WAV::Array{T2,1},waxis::Array{T1,1},w::T1,nz::Int,nx::Int,ext::Int)

# This function computes the source field at the given frequency. The source
# field is vectorized using column major order. This code is intended for
# simultaneous sources.
#
# INPUTS:     isz       - z-index for the source locations
#             isx       - x-index for the source locations
#             ot        - Activation time for the sources
#             WAV       - Source wavelet (frequency domain)
#             waxis     - Angular frequency axis; contains the angular frequency for each entry in 'WAV'
#             w         - Angular frequency to perform modeling for
#             nz        - Number of grid points in the z-direction
#             nx        - Number of grid points in the x-direction
#             ext       - Thickness (number of grid points) of absorbing boundary region
#
# OUTPUTS:    src       - Source field (vector)

    _,iw = findmin(abs(waxis-w))
    src = complex(zeros(nz*nx))
    for isrc = 1:length(isz)
        k = (isx[isrc]+ext-1)*nz + isz[isrc] + ext
        coeff = -WAV[iw]*exp(-im*w*ot[isrc])
        src[k] = coeff
    end

    return src

end

function Source{T1<:AbstractFloat,T2<:Complex}(isz::Int,isx::Int,ot::T1,WAV::Array{T2,1},waxis::Array{T1,1},w::T1,nz::Int,nx::Int,ext::Int)

# This function computes the source field at the given frequency. The source
# field is vectorized using column major order. This code is intended for only
# a single source.
#
# INPUTS:     isz       - z-index for the source location
#             isx       - x-index for the source location
#             ot        - Activation time for the source
#             WAV       - Source wavelet (frequency domain)
#             waxis     - Angular frequency axis; contains the angular frequency for each entry in 'WAV'
#             w         - Angular frequency to perform modeling for
#             nz        - Number of grid points in the z-direction
#             nx        - Number of grid points in the x-direction
#             ext       - Thickness (number of grid points) of absorbing boundary region
#
# OUTPUTS:    src       - Source field (vector)

    _,iw = findmin(abs(waxis-w))
    k = (isx+ext-1)*nz + isz + ext
    coeff = -WAV[iw]*exp(-im*w*ot)

    src = complex(zeros(nz*nx))
    src[k] = coeff

    return src

end