text:
# using Revise
include("D:\\projects\\vsCode\\MedPipe\\MedPipe\\src\\includeAll.jl")

# using Main.LoadFromMonai, Main.HDF5saveUtils,Main.visualizationFromHdf5, Main.distinctColorsSaved
# using PythonCall
# using MedEye3d
# import MedEye3d.ForDisplayStructs.TextureSpec
# using Distributions
# using Clustering
# using IrrationalConstants
# using ParallelStencil
# using HDF5

# #only for finding paths in folder better to do it with pure Julia
# glob=pyimport("glob")
# os=pyimport("os")

# #representing number that is the patient id in this dataset
# patienGroupName="0"
# z=7# how big is the area from which we collect data to construct probability distributions
# klusterNumb = 5# number of clusters - number of probability distributions we will use

# #directory of folder with files in this directory all of the image files should be in subfolder volumes 0-49 and labels labels if one ill use lines below
# data_dir = "D:\\dataSets\\CTORGmini\\"
# #directory where we want to store our HDF5 that we will use
# pathToHDF5="D:\\dataSets\\forMainHDF5\\smallLiverDataSet.hdf5"

# #loading paths of images in our dataset
# train_images = sort(pyconvert(Array,glob.glob(os.path.join(data_dir, "volumes 0-49", "*.nii.gz"))))
# train_labels = sort(pyconvert(Array,glob.glob(os.path.join(data_dir, "labels", "*.nii.gz"))))

# #zipping so we will have tuples with image and label names
# zipped= collect(zip(train_images,train_labels))
# tupl=zipped[1]
# #proper loading
# loaded = LoadFromMonai.loadByMonaiFromImageAndLabelPaths(tupl[1],tupl[2])
# #now We open the hdf5 and save the array with some required metadata

# #!!!!!!!!!! important if you are just creating the hdf5 file  do it with "w" option otherwise do it with "r+"
# fid = h5open(pathToHDF5, "w")
# #fid = h5open(pathToHDF5, "r+")
# gr= getGroupOrCreate(fid, patienGroupName)

# #for this particular example we are intrested only in liver so we will keep only this
# labelArr=map(entry-> UInt8(entry==1),loaded[2])

# #we save loaded and trnsformed data into HDF5 to avoid doing preprocessing every time
# saveMaskBeforeVisualization(fid,patienGroupName,loaded[1],"image", "CT" )
# saveMaskBeforeVisualization(fid,patienGroupName,labelArr,"labelSet", "boolLabel" )

# # here we did default transformations so voxel dimension is set to 1,1,1 in any other case one need to set spacing attribute manually to proper value
# # spacing can be found in metadata dictionary that is third entry in loadByMonaiFromImageAndLabelPaths output
# # here metadata = loaded[3]
# writeGroupAttribute(fid,patienGroupName, "spacing", [1,1,1])

# #******************for display
# #just needed so we will not have 2 same colors for two diffrent informations
# listOfColorUsed= falses(18)

# ##below we define additional arrays that are not present in original data but will be needed for annotations and storing algorithm output

# #manual Modification array
# manualModif = TextureSpec{UInt8}(# choosing number type manually to reduce memory usage
#     name = "manualModif",
#     color =getSomeColor(listOfColorUsed)# automatically choosing some contrasting color
#     ,minAndMaxValue= UInt8.([0,1]) #important to keep the same number type as chosen at the bagining
#     ,isEditable = true ) # we will be able to manually modify this array in a viewer

# algoVisualization = TextureSpec{Float32}(
#     name = "algoOutput",
#     # we point out that we will supply multiple colors
#     isContinuusMask=true,
#     colorSet = [getSomeColor(listOfColorUsed),getSomeColor(listOfColorUsed)]
#     ,minAndMaxValue= Float32.([0,1])# values between 0 and 1 as this represent probabilities
#    )

#     addTextSpecs=Vector{TextureSpec}(undef,2)
#     addTextSpecs[1]=manualModif
#     addTextSpecs[2]=algoVisualization

# #2) primary display of chosen image
# mainScrollDat= loadFromHdf5Prim(fid,patienGroupName,addTextSpecs,listOfColorUsed)

# # now we can use manualModif array to  create annotations

# #****************** constructing probability distributions

# #manually modify in annotator

# ##coordinates of manually set points
# coordsss= GaussianPure.getCoordinatesOfMarkings(eltype(image),eltype(manualModif),  manualModif, image) |>
#     (seedsCoords) ->GaussianPure.getPatchAroundMarks(seedsCoords,z ) |>
#     (patchCoords) ->GaussianPure.allNeededCoord(patchCoords,z )

# #getting patch statistics - mean and covariance
# patchStats = GaussianPure.calculatePatchStatistics(eltype(image),Float64, coordsss, image)

# #separate distribution for each marked point
# distribs = map(patchStat-> fit(MvNormal, reduce(hcat,(patchStat)))  , patchStats  )

# #in order to reduce computational complexity  we will reduce the number of used distributions using kl divergence

# #we are comparing all distributions
# klDivs =map(outerDist->    map(dist->kldivergence( outerDist  ,dist), distribs  ), distribs  )
# klDivsInMatrix = reduce(hcat,(klDivs))
# #clustering with kmeans
# R = kmeans(klDivsInMatrix, klusterNumb; maxiter=200, display=:iter)

# #now identify indexes for some example distributions from each cluster
# indicies = zeros(Int64,klusterNumb )
# a = assignments(R) # get the assignments of points to clusters
# for i in 1:klusterNumb
#     for j in 1:length(distribs)
#         if(a[j] == i)
#             indicies[i]=j
#         end
#     end
# end
# indicies

# #ditributions from diffrent clusters
# chosenDistribs = map(ind->distribs[ind] ,indicies)

# # imageDataset.attrs['dataType']=  "CT"
# # labelBoolDataset.attrs['dataType']=  "boolLabel"

# # using Revise,CUDA

# # using MedEye3d
# # import MedEye3d.ForDisplayStructs.TextureSpec
# # using Distributions
# # using Clustering
# # using ParallelStencil

# # import MedEye3d
# # import MedEye3d.ForDisplayStructs
# # import MedEye3d.ForDisplayStructs.TextureSpec
# # using ColorTypes
# # import MedEye3d.SegmentationDisplay

# # import MedEye3d.DataStructs.ThreeDimRawDat
# # import MedEye3d.DataStructs.DataToScrollDims
# # import MedEye3d.DataStructs.FullScrollableDat
# # import MedEye3d.ForDisplayStructs.KeyboardStruct
# # import MedEye3d.ForDisplayStructs.MouseStruct
# # import MedEye3d.ForDisplayStructs.ActorWithOpenGlObjects
# # import MedEye3d.OpenGLDisplayUtils
# # import MedEye3d.DisplayWords.textLinesFromStrings
# # import MedEye3d.StructsManag.getThreeDims

# # listOfColorUsed= falses(18)

# # CTIm= loaded[1]
# # labell=loaded[2]

# # #manual Modification array
# # CTImm= TextureSpec{Float32}(
# #     name= "CTIm",
# #     numb= Int32(3),
# #     isMainImage = true,
# #     minAndMaxValue= Float32.([0,100]))

# # algoVisualization = TextureSpec{Float32}(
# #     name = "algoOutput",
# #     # we point out that we will supply multiple colors
# #     isContinuusMask=true,
# #     colorSet = [getSomeColor(listOfColorUsed),getSomeColor(listOfColorUsed)]
# #     ,minAndMaxValue= Float32.([0,1])# values between 0 and 1 as this represent probabilities
# #    )

# #    manualModif=  TextureSpec{UInt8}(
# #     name = "manualModif",
# #     numb= Int32(2),
# #     color = RGB(0.0,1.0,0.0)
# #     ,minAndMaxValue= UInt8.([0,1])
# #     ,isEditable = true
# #    )

# #     addTextSpecs=Vector{TextureSpec}(undef,3)
# #     addTextSpecs[1]=manualModif
# #     addTextSpecs[2]=CTImm
# #     addTextSpecs[3]=algoVisualization

# #     datToScrollDimsB= MedEye3d.ForDisplayStructs.DataToScrollDims(imageSize=  size(CTIm) ,voxelSize=(1,1,1), dimensionToScroll = 3 );

# #     import MedEye3d.DisplayWords.textLinesFromStrings

# #     mainLines= textLinesFromStrings(["main Line1", "main Line 2"]);
# #     supplLines=map(x->  textLinesFromStrings(["sub  Line 1 in $(x)", "sub  Line 2 in $(x)"]), 1:size(CTIm)[3] );

# #     import MedEye3d.StructsManag.getThreeDims

# #     tupleVect = [("algoOutput",labell) ,("CTIm",CTIm),("manualModif",zeros(UInt8,size(CTIm)) ) ]
# #     slicesDat= getThreeDims(tupleVect )

# #     fractionOfMainIm= Float32(0.8);

# #     mainScrollDat = FullScrollableDat(dataToScrollDims =datToScrollDimsB
# #     ,dimensionToScroll=1 # what is the dimension of plane we will look into at the beginning for example transverse, coronal ...
# #     ,dataToScroll= slicesDat
# #     ,mainTextToDisp= mainLines
# #     ,sliceTextToDisp=supplLines );

# #     SegmentationDisplay.coordinateDisplay(addTextSpecs ,fractionOfMainIm ,datToScrollDimsB ,1000);

# #     Main.SegmentationDisplay.passDataForScrolling(mainScrollDat);


text:
# Copyright (C) 2020 Argonne National Laboratory
# Written by Alinson Santos Xavier <axavier@anl.gov>

using RELOG, Cbc, JuMP, Printf, JSON, MathOptInterface.FileFormats

@testset "build" begin
    basedir = dirname(@__FILE__)
    instance = RELOG.parsefile("$basedir/../../instances/s1.json")
    graph = RELOG.build_graph(instance)
    model = RELOG.build_model(instance, graph, Cbc.Optimizer)
    set_optimizer_attribute(model, "logLevel", 0)

    process_node_by_location_name =
        Dict(n.location.location_name => n for n in graph.process_nodes)

    shipping_node_by_loc_and_prod_names = Dict(
        (n.location.location_name, n.product.name) => n for n in graph.plant_shipping_nodes
    )

    @test length(model[:flow]) == 76
    @test length(model[:dispose]) == 16
    @test length(model[:open_plant]) == 12
    @test length(model[:capacity]) == 12
    @test length(model[:expansion]) == 12

    l1 = process_node_by_location_name["L1"]
    v = model[:capacity][l1, 1]
    @test lower_bound(v) == 0.0
    @test upper_bound(v) == 1000.0

    v = model[:expansion][l1, 1]
    @test lower_bound(v) == 0.0
    @test upper_bound(v) == 750.0

    v = model[:dispose][shipping_node_by_loc_and_prod_names["L1", "P2"], 1]
    @test lower_bound(v) == 0.0
    @test upper_bound(v) == 1.0
end

ext: jl

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
"""
    module Geometry

This module defines basic geometrical objects.
"""
module Geometry

import Base: ==, length, in, iterate, getindex, lastindex, isempty, eltype, keys
import Base: union, setdiff, intersect, issubset

using DocStringExtensions
using StaticArrays
using LinearAlgebra
using ForwardDiff # for computing derivatives of parametric elements
using OrderedCollections

using WaveProp.Utils

include("point.jl")
include("hyperrectangle.jl")
include("referenceshapes.jl")
include("entities.jl")
include("domain.jl")
include("parametricbody.jl")

export
    # abstract types
    AbstractReferenceShape,
    AbstractEntity,
    AbstractParametricBody,
    # types
    ElementaryEntity,
    ParametricEntity,
    ParametricCurve,
    ClosedEntity,
    Circle,
    Kite,
    Domain,
    ReferenceLine,
    ReferenceTriangle,
    ReferenceTetrahedron,
    ReferenceSquare,
    HyperRectangle,
    SVector,
    # functions
    entities,
    key,
    remove,
    assertequaldim,
    boundary,
    skeleton,
    internal_boundary,
    external_boundary,
    domain,
    parametrization,
    jacobian,
    normal,
    ambient_dimension,
    geometric_dimension,
    bounding_box,
    center,
    diameter,
    radius,
    low_corner,
    high_corner,
    measure,
    number_of_nodes,
    vertices,
    line,
    # global variables
    TAGS,
    ENTITIES

end # module

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module XLParser
    export @xlexpr_str
    export @xl_str
    export xlexpr

    using Tokenize
    using Espresso

    # we bait and switch string concatenation (ampersand in xl) with the pipe
    # character, since it is placed on the appropriate precedence level in Julia:
    # https://github.com/JuliaLang/julia/blob/master/src/julia-parser.scm
    dummy_ampersand = :(|>)

    function _sanitize_excel_formula_strings(str)
        chars = collect(str)
        str_builder = Array{Union{String, Char}, 1}()

        i = 0
        inquote = false
        while (i+=1) <= length(chars)

            # Mark the beginning of a string
            if !inquote && chars[i] == '"'
                inquote = true
                push!(str_builder, "raw\"")

            # escaping backslashes if they are before a quote (double them up)
            elseif inquote && chars[i] == '\\'
                cnt = 1
                while (i+=1) <= length(chars) && chars[i] == '\\'
                    cnt+=1
                end

                if chars[i] == '"'
                    push!(str_builder, "\\"^(cnt*2))
                else
                    push!(str_builder, "\\"^(cnt))
                end

                # reset our peek-forward
                i -= 1

            # escaping of a quote
            elseif inquote && (i+1) <= length(chars) && chars[i] == '"' && chars[i+1] == '"'
                push!(str_builder, "\\\"")
                i+=1

            # closing the string
            elseif inquote && chars[i] == '"'
                push!(str_builder, '"')
                inquote = false

            # normal passthrough
            else
                push!(str_builder, chars[i])
            end
        end

        return join(str_builder, "")
    end

    function _sanitize_excel_formula_tokens(str)
        # ensure excel strings are converted to Julia syntax
        str = _sanitize_excel_formula_strings(str)
        str = replace(str, "\n" => " ")

        tokens = collect(tokenize(str))

        i = 0
        while (i+=1) <= length(tokens)

            # convert names to lowercase
            if tokens[i].kind == Tokenize.Tokens.IDENTIFIER && tokens[i].val != lowercase(tokens[i].val)
                tokens[i] = first(tokenize(lowercase(tokens[i].val)))
                i-=1

            # convert if KEYWORD token to xl_if function token
            elseif tokens[i].kind == Tokenize.Tokens.IF
                tokens[i] = first(tokenize("xl_if"))

            # convert <> to !=
            elseif (i != length(tokens) &&
                tokens[i].kind == Tokenize.Tokens.LESS &&
                tokens[i+1].kind == Tokenize.Tokens.GREATER)

                tokens[i] = first(tokenize("!="))
                tokens[i+1] = first(tokenize(""))

            # convert = to ==
            elseif (tokens[i].kind == Tokenize.Tokens.EQ)
                tokens[i] = first(tokenize("=="))

            # add whitespace between repeating minus signs to avoid syntax error
            elseif (tokens[i].token_error == Tokenize.Tokens.INVALID_OPERATOR &&
                        tokens[i].val == "--")
                tokens[i] = first(tokenize("-"))
                insert!(tokens, i+1, first(tokenize(" ")))
                insert!(tokens, i+1, first(tokenize("-")))
                insert!(tokens, i+1, first(tokenize(" ")))

            # fill empty commas: FOO(,,) -> FOO(nothing,nothing,nothing)
            elseif (tokens[i].kind == Tokenize.Tokens.COMMA)
                # test if left of comma is bracket
                j = i-1
                if j >= 1 && tokens[j].kind == Tokenize.Tokens.WHITESPACE
                    j -= 1
                end

                if j >= 1 && tokens[j].kind == Tokenize.Tokens.LPAREN
                    insert!(tokens, j+1, first(tokenize("nothing")))
                    i += 1
                end

                # test if right of comma is another comma or )
                j = i+1
                if j <= length(tokens) && tokens[j].kind == Tokenize.Tokens.WHITESPACE
                    j += 1
                end

                if j <= length(tokens) && (tokens[j].kind == Tokenize.Tokens.COMMA ||
                                           tokens[j].kind == Tokenize.Tokens.RPAREN)
                    insert!(tokens, j, first(tokenize("nothing")))
                    i += 1
                end

            # Replace ampersand with a lower precedence infix operator
            elseif  tokens[i].kind == Tokenize.Tokens.AND
                tokens[i] = first(tokenize(string(dummy_ampersand)))
            end
        end

        return join(untokenize.(tokens), "")
    end

    function xlexpr(str)
        # turn string into a valid (but logically wrong) Julia expression
        parsable_string = _sanitize_excel_formula_tokens(str)
        xl = Meta.parse(parsable_string)

        # now correct this logically wrong expression into a correct one

        # index syntax, eg: a[1:5]
        xl = rewrite_all(xl, :(index(_a,_b):index(_a,_d)), :(_a[_b:_d]))
        xl = rewrite_all(xl, :(index(_a,_b)), :(_a[_b]))

        #TODO: test expr contains index(_a,_b):index(_c,_d) and throw an error if so
        xl = rewrite_all(xl, :(index(_a,_b):index(_a,_d)), :(_a[_b:_d]))

        # replace & infix with string function
        xl = rewrite_all(xl, dummy_ampersand, :(string))
        xl = rewrite_all(xl, :(string(string(_a...), _b)),
                             :(string(_a..., _b)))

        # replace and(1,2,3) with 1 && (2 && 3)
        while xl != (xl_new = rewrite_all(xl, :(and(_a, _b...)),
                                              :(_a && and(_b...))))
            xl = xl_new
        end
        xl = rewrite_all(xl, :(and(_a)), :(_a))

        # replace or(1,2,3) with 1 || (2 || 3)
        while xl != (xl_new = rewrite_all(xl, :(or(_a, _b...)),
                                              :(_a || or(_b...))))
            xl = xl_new
        end
        xl = rewrite_all(xl, :(or(_a)), :(_a))

        # replace not(1) with !(1)
        xl = rewrite_all(xl, :(not(_a)), :(!(_a)))

        # replace excel if functions with julia if statements
        xl = rewrite_all(xl, :(xl_if(_a,_b,_c)), :((_a ? _b : _c)))
        xl = rewrite_all(xl, :(xl_if(_a,_b)), :((_a ? _b : false)))

        return xl
    end

    macro xlexpr_str(str)
        # quot allows macro to not execute but to return the expression
        Meta.quot(xlexpr(str))
    end

    macro xl_str(str)
        # esc allows macro to execute in caller's scope
        esc(xlexpr(str))
    end

end


text:
"""
Aligns a collection of functions using the elastic square-root slope (srsf)
framework.

    srsf_align(f, timet; method="mean", smooth=false, sparam=10, lam=0.0,
               optim="DP", MaxItr=20)
    :param f: array of shape (M,N) of N functions with M samples
    :param timet: vector of size M describing the sample points
    :param method: (string) calculate Karcher Mean or Median
    (options = "mean" or "median") (default="mean")
    :param smooth: Smooth the data using a box filter (default = false)
    :param sparam: Number of times to run smoothing filter (default 10)
    :param lam: controls the elasticity (default = 0)
    :param optim: optimization method to find warping, default is
                  Simultaneous Alignment ("SIMUL"). Other options are
                  Dynamic Programming ("DP2"), Riemannian BFGS
                  ("RBFGS")
    :param MaxItr: Maximum number of iterations

    Returns fdawarp type containing
    :return fn: aligned functions - array of shape (M,N) of N
                functions with M samples
    :return qn: aligned srsfs - similar structure to fn
    :return q0: original srsfs - similar structure to fn
    :return fmean: function mean or median - vector of length N
    :return mqn: srvf mean or median - vector of length N
    :return gam: warping functions - similar structure to fn
    :return orig_var: Original Variance of Functions
    :return amp_var: Amplitude Variance
    :return phase_var: Phase Variance
"""
function srsf_align(f, timet; method="mean", smooth=false, sparam=10, lam=0.0,
                    optim="DP", MaxItr=20)
    M, N = size(f);
    if smooth
        smooth_data!(f, sparam);
    end

    if M > 500
        parallel = true;
    elseif N > 100
        parallel = true;
    else
        parallel = false;
    end

    epsilon = eps(Float64);
    f0 = copy(f);

    methods_ = ["mean", "median"];

    method_ = "";
    for ii in methods_
        if method == ii
            method_ = ii;
        end
    end
    if method_ == ""
        warn("improper method, defaulting to mean")
        method = "mean";
    end

    # Compute SRSF
    binsize = mean(diff(timet));
    q = Array{Float64}(undef, (M, N));
    for ii in 1:N
        q[:, ii] = f_to_srsf(f[:, ii], timet);
    end
    fo = vec(f[1,:]);
    println("Initializing...")
    mnq = mean(q,dims=2);
    d1 = repeat(mnq,1,N);
    d = (q - d1).^2;
    dqq = sqrt.(sum(d,dims=1));
    min_ind = argmax(vec(dqq));
    mq = q[:, min_ind];
    mf = f[:, min_ind];

    if parallel
        gam = @distributed (hcat) for i=1:N
            optimum_reparam(mq, timet, q[:, i], lam, method=optim,
                            f1o=mf[1], f2o=fo[i]);
        end
    else
        gam = optimum_reparam(mq, timet, q, lam, method=optim,
                              f1o=mf[1], f2o=fo);
    end

    gamI = sqrt_mean_inverse(gam);
    xout = (timet[end] -timet[1]) .* gamI .+ timet[1];
    mf = approx(timet, mf, xout);
    mq = f_to_srsf(mf ,timet);

    # Compute Karcher Mean
    if method == "mean"
        @printf("Compute Karcher Mean of %d functions in SRSF space..\n",N)
    elseif method == "median"
        @printf("Compute Karcher Median of %d functions in SRSF space..\n",N)
    end

    ds = zeros(MaxItr+2);
    ds[1] = Inf;
    qun = zeros(MaxItr+1);
    tmp = zeros(M, MaxItr+2);
    tmp[:, 1] = mq;
    mq = copy(tmp);
    tmp = zeros(M, MaxItr+2);
    tmp[:, 1] = mf;
    mf = copy(tmp);
    tmp = zeros(M, N, MaxItr+2);
    tmp[:, :, 1] = f;
    f = copy(tmp);
    tmp = zeros(M, N, MaxItr+2);
    tmp[:, :, 1] = q;
    q = copy(tmp);
    r1 = 0;

    for r in 1:MaxItr
        @printf("updating step: r=%d\n", r)
        if r == MaxItr
            println("maximal number of iterations reached")
        end

        # Matching Step
        if parallel
            gam = @distributed (hcat) for i=1:N
                optimum_reparam(mq[:,r], timet, q[:, i, 1], lam, method=optim,
                                f1o=mf[1,r], f2o=fo[i]);
            end
        else
            gam = optimum_reparam(mq[:,r], timet, q[:,:,1], lam, method=optim,
                                  f1o=mf[1,r], f2o=fo);
        end

        gam_dev = zeros(M,N);
        for k in 1:N
            xout = (timet[end]-timet[1]) .* gam[:, k] .+ timet[1];
            f[:, k, r+1] = approx(timet, f[:, k, 1], xout);
            q[:, k, r+1] = f_to_srsf(f[:, k, r+1], timet);
            gam_dev = gradient(gam[:, k], 1/(M-1));
        end

        mqt = mq[:, r];
        a = repeat(mqt,1,N);
        d = (q[:, :, r+1] - a).^2;
        if method == "mean"
            d1 = sum(trapz(timet, d));
            d2 = sum(trapz(timet, (1 .-sqrt.(gam_dev)).^2));
            ds_tmp = d1 + lam * d2;
            ds[r+1] = ds_tmp;

            # Minimization Step
            # compute the mean of the matched function
            qtemp = q[:, :, r+1];
            mq[:, r+1] = mean(qtemp, dims=2);
            mf[:, r+1] = mean(f[:, :, r+1], dims=2);

            qun[r] = norm(mq[:,r+1] - mq[:,r])/norm(mq[:,r]);
        end

        if method == "median"
            d1 = sqrt(sum(trapz(timet, d)));
            d2 = sum(trapz(timet, (1 .-sqrt.(gam_dev)).^2));
            ds_tmp = d1 * lam * d2;
            ds[r+1] = ds_tmp;

            # Minimization Step
            # compute the median of the matched functions
            dist_iinv = ds[r+1].^(-1);
            qtemp = q[:, :, r+1]/ds[r+1];
            mq[:, r+1] = sum(qtemp,2) .* dist_iinv;
            mf[:, r+1] = sum(f[:,:,r+1]/ds[r+1],2) .* dist_iinv;

            qun[r] = norm(mq[:, r+1] - mq[:,r])/ norm(mq[:,r]);
        end

        r1 = r;
        if (qun[r] < 1e-2) || (r >= MaxItr)
            break
        end
    end

    # Last Step with Centering of gam
    r = r1 + 1;
    if parallel
        gam = @distributed (hcat) for i=1:N
            optimum_reparam(mq[:,r], timet, q[:, i, 1], lam, method=optim,
                            f1o=mf[1,r], f2o=fo[i]);
        end
    else
        gam = optimum_reparam(mq[:,r], timet, q[:,:,1], lam, method=optim,
                              f1o=mf[1,r], f2o=fo);
    end

    gam_dev = zeros(M,N);
    for k in 1:N
        gam_dev = gradient(gam[:, k], 1/(M-1));
    end

    gamI = sqrt_mean_inverse(gam);
    gamI_dev = gradient(gamI, 1/(M-1));
    timet0 = (timet[end]-timet[1]) .* gamI  .+ timet[1];
    mq[:, r+1] = approx(timet, mq[:, r], timet0) .* sqrt.(gamI_dev);

    for k in 1:N
        q[:, k, r+1] = approx(timet, q[:, k, r+1], timet0) .* sqrt.(gamI_dev);
        f[:, k, r+1] = approx(timet, f[:, k, r+1], timet0);
        gam[:, k] = approx(timet, gam[:, k], timet0);
    end

    # Aligned data & stats
    fn = f[:, :, r];
    qn = q[:, :, r];
    q0 = q[:, :, 1];
    mean_f0 = mean(f0, dims=2);
    std_f0 = std(f0, dims=2);
    mean_fn = mean(fn, dims=2);
    std_fn = std(fn, dims=2);
    mqn = mq[:, r];
    fmean = mean(f0[1,:]) .+ cumtrapz(timet, mqn .* abs.(mqn));

    fgam = zeros(M, N);
    for ii in 1:N
        xout = (timet[end] -timet[1]) .* gam[:, ii] .+ timet[1];
        fgam[:, ii] = approx(timet, fmean, xout);
    end
    var_fgam = var(fgam,dims=2);

    orig_var = trapz(timet, std_f0.^2);
    amp_var = trapz(timet, std_fn.^2);
    phase_var = trapz(timet, var_fgam);

    out = fdawarp(f0, timet, fn, qn, q0, fmean, mqn, gam,
                  orig_var, amp_var, phase_var, qun, lam,
                  method, optim, gamI, false, Array{Float64}(undef,0),
                  Array{Float64}(undef,0), Array{Float64}(undef,0),
                  Array{Float64}(undef,0))
    return out
end

"""
Aligns a collection of functions while extracting principal components.
The functions are aligned to the principal components

    align_fPCA(f, timet; num_comp=3, smooth=false, sparam=10, MaxItr=50)
    :param f: array of shape (M,N) of N functions with M samples
    :param timet: vector of size M describing the sample points
    :param num_comp: Number of components (default = 3)
    :param smooth: Smooth the data using a box filter (default = false)
    :param sparam: Number of times to run smoothing filter (default 10)
    :param MaxItr: Maximum number of iterations

    Returns Dict containing
    :return fn: aligned functions - array of shape (M,N) of N functions with M
                samples
    :return qn: aligned srvfs - similar structure to fn
    :return q0: original srvf - similar structure to fn
    :return mqn: srvf mean or median - vector of length M
    :return gam: warping functions - similar structure to fn
    :return q_pca: srsf principal directions
    :return f_pca: functional principal directions
    :return latent: latent values
    :return coef: coefficients
    :return U: eigenvectors
    :return orig_var: Original Variance of Functions
    :return amp_var: Amplitude Variance
    :return phase_var: Phase Variance
"""
function align_fPCA(f, timet; num_comp=3, smooth=false, sparam=10, MaxItr=50)
    lam = 0.0;
    coef = collect(-2:2);
    Nstd = length(coef);
    M, N = size(f);

    if smooth
        smooth_data!(f, sparam);
    end

    if M > 500
        parallel = true;
    elseif N > 100
        parallel = true;
    else
        parallel = false;
    end

    epsilon = eps(Float64);
    f0 = copy(f);

    # Compute SRSF
    binsize = mean(diff(timet));
    q = Array{Float64}(undef, (M, N));
    for ii in 1:N
        q[:, ii] = f_to_srsf(f[:, ii], timet);
    end
    println("Initializing...")
    mnq = mean(q,dims=2);
    d1 = repeat(mnq,1,N);
    d = (q - d1).^2;
    dqq = sqrt.(sum(d,dims=1));
    min_ind = argmax(vec(dqq));
    @printf("Compute %d functions in SRSF space to %d fPCA components..\n",N,num_comp)

    # Compute Karcher Mean
    itr = 1;
    mq = zeros(M, MaxItr+1);
    mq[:, 1] = q[:, min_ind];
    fi = zeros(M, N, MaxItr+1);
    fi[:, :, 1] = f;
    qi = zeros(M, N, MaxItr+1);
    qi[:, :, 1] = q;
    gam = zeros(M, N, MaxItr+1);
    cost = zeros(MaxItr+1);
    itrf = 0;

    while itr<=MaxItr
        @printf("updating step: r=%d\n", itr)
        if itr == MaxItr
            println("maximal number of iterations reached")
        end

        # PCA Step
        a = repeat(mq[:, itr], 1, N);
        qhat_cent = qi[:, :, itr] - a;
        xmean = vec(mean(qi[:, :, itr],dims=2));
        K = Statistics.covm(qi[:, :, itr], xmean, 2);
        U, s, V = svd(K);

        alpha_i = zeros(num_comp, N);
        for ii in 1:num_comp
            for jj in 1:N
                alpha_i[ii,jj] = trapz(timet, qhat_cent[:, jj] .* U[:, ii]);
            end
        end

        U1 = U[:,1:num_comp];
        tmp = U1 * alpha_i;
        qhat = a + tmp;

        # Matching Step
        if parallel
            gam_t = @distributed (hcat) for i=1:N
                optimum_reparam(qhat[:,i], timet, qi[:,i,itr], lam,
                                method="DP");
            end
            gam[:, :, itr] = gam_t;
        else
            gam[:, :, itr] = optimum_reparam(qhat, timet, qi[:,:,itr], lam,
                                             method="DP");
        end

        for k in 1:N
            xout = (timet[end] -timet[1]) .* gam[:, k, itr] .+ timet[1];
            fi[:, k, itr+1] = approx(timet, fi[:, k, itr], xout);
            qi[:, k, itr+1] = f_to_srsf(fi[:, k, itr+1], timet);
        end

        qtemp = qi[:, :, itr+1];
        mq[:, itr+1] = mean(qtemp, dims=2);

        cost_temp = zeros(N);
        for ii in 1:N
            cost_temp[ii] = norm(qtemp[:,ii] - qhat[:,ii]).^2;
        end

        cost[itr+1] = mean(cost_temp);

        if abs(cost[itr+1] - cost[itr] < 1e-6)
            break
        end

        itr += 1;
        itrf = itr;
    end

    if itrf >= MaxItr
        itrf = MaxItr;
    else
        itrf += 1;
    end
    cost = cost[1:itrf];

    # Aligned data & stats
    fn = fi[:, :, itrf];
    qn = qi[:, :, itrf];
    q0 = qi[:, :, 1];
    mean_f0 = mean(f0, dims=2);
    std_f0 = std(f0, dims=2);
    mqn = mq[:, itrf];
    gamf = gam[:, :, 1];
    for k in 1:itrf-1
        gam_k = gam[:, :, k];
        for l in 1:N
            time0 = (timet[end]-timet[1]) .* gam_k[:, l] .+ timet[1];
            gamf[:,l] = approx(timet, gamf[:,l], time0);
        end
    end
    fmean = mean(f0[1,:]) .+ cumtrapz(timet, mqn .* abs.(mqn));

    # Center Mean
    gamI = sqrt_mean_inverse(gamf);
    gamI_dev = gradient(gamI, 1/(M-1));
    timet0 = (timet[end] -timet[1]) .* gamI  .+ timet[1];
    mqn = approx(timet, mqn, timet0) .* sqrt.(gamI_dev);
    for k in 1:N
        qn[:, k] = approx(timet, qn[:, k], timet0) .* sqrt.(gamI_dev);
        fn[:, k] = approx(timet, fn[:, k], timet0);
        gamf[:, k] = approx(timet, gamf[:, k], timet0);
    end
    mean_fn = mean(fn, dims=2);
    std_fn = std(fn, dims=2);

    fgam = zeros(M, N);
    for ii in 1:N
        xout = (timet[end] - timet[1]) .* gamf[:, ii] .+ timet[1];
        fgam[:, ii] = approx(timet, fmean, xout);
    end

    var_fgam = var(fgam,dims=2);
    orig_var = trapz(timet, std_f0.^2);
    amp_var = trapz(timet, std_fn.^2);
    phase_var = trapz(timet, var_fgam);

    # Get Final PCA
    out_warp = fdawarp(f0, timet, fn, qn, q0, fmean, mqn, gamf,
                       orig_var, amp_var, phase_var, cost, lam,
                       "PCA", "DP", gamI, false, Array{Float64}(undef,0),
                       Array{Float64}(undef,0), Array{Float64}(undef,0),
                       Array{Float64}(undef,0))

    out = vert_fPCA(out_warp, no=num_comp);

    out2 = Dict("fn" => fn, "qn" => qn, "q0" => q0, "fmean" => fmean,
                "mqn" => mqn, "gam" => gamf, "q_pca" => out.q_pca,
                "f_pca" => out.f_pca, "latent" => out.latent,
                "coef" => out.coef, "U" => out.U, "orig_var" => orig_var,
                "amp_var" => amp_var, "phase_var" => phase_var);

    return out2
end


text:
module WavesAndEigenvalues

include("Meshutils.jl")
include("NLEVP.jl")
include("Helmholtz.jl")
include("APE.jl")
include("network.jl")

#using .Meshutils
#include("meshutils_export.jl")

end # module


text:
"""
    RandomSelector(seed, holdout_fraction=1//2, block_size=Day(1), offset=Day(0))

Determine holdout set by randomly subsampling contiguous blocks of size `block_size`
without replacement using a `MersenneTwister` seeded with `seed`.
The probability of any given block being in the holdout set is given by `holdout_fraction`.

The `offset` is rarely needed, but is used to control block boundries.
It is given as a offset relative to _Monday 1st Jan 1900_.
For example, with the default offset of `Day(0)`, and if using a `Week(1)` `block_size`,
then every block will start on a Monday, and will go for 1 or more weeks from there.

Note that at the boundries of the partitioned dates the blocks may not be of size
`block_size` if they go over the edge -- this is infact the common case.
"""
struct RandomSelector <: DateSelector
    seed::Int
    holdout_fraction::Real
    block_size::DatePeriod
    offset::DatePeriod

    function RandomSelector(seed, holdout_fraction=1//2, block_size=Day(1), offset=Day(0))
        if !(0 <= holdout_fraction <= 1)
            throw(DomainError(
                holdout_fraction,
                "holdout fraction must be between 0 and 1 (inclusive)"
            ))
        end
        if block_size < Day(1)
            throw(DomainError(block_size, "block_size must be at least 1 day."))
        end
        return new(seed, holdout_fraction, block_size, offset)
    end
end

function Iterators.partition(dates::StepRange{Date, Day}, s::RandomSelector)
    sd, ed = extrema(dates)

    rng = MersenneTwister(s.seed)

    holdout_dates = Date[]
    initial_time = _initial_date(s, dates)
    curr_window = initial_time:step(dates):(initial_time + s.block_size - step(dates))
    while first(curr_window) <= ed
        # Important: we must generate a random number for every block even before the start
        # so that the `rng` state is updated constistently no matter when the start is
        # and thus `partition` is invarient on the start date
        r = rand(rng)

        # optimization: only creating holdout window if intersect not empty
        if last(curr_window) >= sd
            if r < s.holdout_fraction
                curr_active_window = curr_window ∩ dates  # handle being near boundries
                append!(holdout_dates, curr_active_window)
            end
        end
        curr_window = curr_window .+ s.block_size
    end

    return _getdatesets(dates, holdout_dates)
end


text:
#=

Filename: asset_pricing.jl

@authors: Spencer Lyon, Tom Sargent, John Stachurski

Computes asset prices with a Lucas style discount factor when the endowment
obeys geometric growth driven by a finite state Markov chain.  That is,

.. math::
    d_{t+1} = g(X_{t+1}) d_t

where

    * :math:`\{X_t\}` is a finite Markov chain with transition matrix P.

    * :math:`g` is a given positive-valued function

References
----------

    http://quant-econ.net/py/markov_asset.html

=#

using QuantEcon

# A default Markov chain for the state process
rho = 0.9
sigma = 0.02
n = 25
default_mc = tauchen(n, rho, sigma)

type AssetPriceModel
    beta :: Float64    # Discount factor
    gamma :: Float64   # Coefficient of risk aversion
    mc :: MarkovChain  # State process
    n :: Int           # Number of states
    g :: Function      # Function mapping states into growth rates
end

function AssetPriceModel(;beta=0.96, gamma=2.0, mc=default_mc, g=exp)
    n = size(mc.p)[1]
    return AssetPriceModel(beta, gamma, mc, n, g)
end

"""
Stability test for a given matrix Q.
"""
function test_stability(ap::AssetPriceModel, Q::Matrix)
    sr = maximum(abs(eigvals(Q)))
    if sr >= 1 / ap.beta
        msg = "Spectral radius condition failed with radius = $sr"
        throw(ArgumentError(msg))
    end
end

"""
Computes the price-dividend ratio of the Lucas tree.

"""
function tree_price(ap::AssetPriceModel)
    # == Simplify names, set up matrices  == #
    beta, gamma, P, y = ap.beta, ap.gamma, ap.mc.p, ap.mc.state_values
    y = reshape(y, 1, ap.n)
    J = P .* ap.g(y).^(1 - gamma)

    # == Make sure that a unique solution exists == #
    test_stability(ap, J)

    # == Compute v == #
    I = eye(ap.n)
    Ones = ones(ap.n)
    v = (I - beta * J) \ (beta * J * Ones)

    return v
end

"""
Computes price of a consol bond with payoff zeta

"""
function consol_price(ap::AssetPriceModel, zeta::Float64)
    # == Simplify names, set up matrices  == #
    beta, gamma, P, y = ap.beta, ap.gamma, ap.mc.p, ap.mc.state_values
    y = reshape(y, 1, ap.n)
    M = P .* ap.g(y).^(- gamma)

    # == Make sure that a unique solution exists == #
    test_stability(ap, M)

    # == Compute price == #
    I = eye(ap.n)
    Ones = ones(ap.n)
    p = (I - beta * M) \ ( beta * zeta * M * Ones)

    return p
end

"""
Computes price of a perpetual call option on a consol bond.

"""
function call_option(ap::AssetPriceModel, zeta::Float64, p_s::Float64, epsilon=1e-7)

    # == Simplify names, set up matrices  == #
    beta, gamma, P, y = ap.beta, ap.gamma, ap.mc.p, ap.mc.state_values
    y = reshape(y, 1, ap.n)
    M = P .* ap.g(y).^(- gamma)

    # == Make sure that a unique console price exists == #
    test_stability(ap, M)

    # == Compute option price == #
    p = consol_price(ap, zeta)
    w = zeros(ap.n, 1)
    error = epsilon + 1
    while (error > epsilon)
        # == Maximize across columns == #
        w_new = max(beta * M * w, p - p_s)
        # == Find maximal difference of each component and update == #
        error = maximum(abs(w-w_new))
        w = w_new
    end

    return w
end


text:
using Wilsonloop
import Wilsonloop:DwDU,get_direction,get_position#,Adjoint_GLink
import Gaugefields:evaluate_gaugelinks!,clear_U!,set_wing_U!

#=
g1,g2,g3 = mk_gamma(1)
const γmatrices = g1
const γ5 = γmatrices[:,:,5]

const rplusγ1 = g2
const rminusγ1 = g3
=#

function shift_wilsonloop(w::Wilsonline{Dim},shift) where Dim
    wshifted = Wilsonline(Dim=Dim)
    #show(w)
    for i=1:length(w)
        link = w[i]
        position = collect(get_position(link)) .- shift
        #println("link: ", get_position(link))
        #println("position, ",position)
        link_shifted = Wilsonloop.set_position(link,Tuple(position))
        push!(wshifted,link_shifted)
    end
    return wshifted
end

struct Linkinfo{Dim}
    positions::Vector{NTuple{Dim,Int64}}
    directions::Vector{Int8}
    isdagvectors::Vector{Bool}

    function Linkinfo(line_for_gaugefields::Wilsonline{Dim}) where Dim
        numlinks = length(line_for_gaugefields)
        positions = Array{NTuple{Dim,Int64},1}(undef,numlinks)
        directions = Array{Int8,1}(undef,numlinks)
        isdagvectors = Array{Bool,1}(undef,numlinks)
        for i=1:numlinks
            glink = line_for_gaugefields[i]
            positions[i] = get_position(glink)
            directions[i] = get_direction(glink)
            isdagvectors[i] = isdag(glink)
            #isdagvectors[i] = ifelse(typeof(glink ) <: Adjoint_GLink,true,false)
        end
        return new{Dim}(positions,directions,isdagvectors)
    end
end

struct General_Dirac_term{Dim}
    coefficient::Float64
    matrix_for_spinors::Matrix{ComplexF64}
    matrix_for_spinors_dagger::Matrix{ComplexF64}
    matrix_for_spinors_γ5M::Matrix{ComplexF64}
    line_for_gaugefields::Wilsonline{Dim}
    linkinfo::Linkinfo{Dim}
    #positions::Vector{NTuple{Dim,Int64}}
    #directions::Vector{Int8}
    #isdagvectors::Vector{Bool}
    derivatives::Vector{Vector{DwDU{Dim}}}
    linkinfo_derivatives_right::Vector{Vector{Linkinfo{Dim}}}
    linkinfo_derivatives_leftdag::Vector{Vector{Linkinfo{Dim}}}
    derivatives_dag::Vector{Vector{DwDU{Dim}}}
    linkinfo_derivatives_dag_right::Vector{Vector{Linkinfo{Dim}}}
    linkinfo_derivatives_dag_leftdag::Vector{Vector{Linkinfo{Dim}}}
    startingpoint::NTuple{Dim,Int64}
    endingpoint::NTuple{Dim,Int64}

    function General_Dirac_term(coefficient,matrix_for_spinors,line_for_gaugefields::Wilsonline{Dim},
        startingpoint,
        endingpoint
        ) where Dim

        matrix_for_spinors_dagger = γ5*matrix_for_spinors*γ5
        matrix_for_spinors_γ5M = γ5*matrix_for_spinors

        linkinfo = Linkinfo(line_for_gaugefields)

        derivatives = Array{Vector{DwDU{Dim}},1}(undef,Dim)
        linkinfo_derivatives_right = Vector{Vector{Linkinfo{Dim}}}(undef,Dim)
        linkinfo_derivatives_leftdag = Vector{Vector{Linkinfo{Dim}}}(undef,Dim)

        derivatives_dag = Array{Vector{DwDU{Dim}},1}(undef,Dim)
        linkinfo_derivatives_dag_right = Vector{Vector{Linkinfo{Dim}}}(undef,Dim)
        linkinfo_derivatives_dag_leftdag = Vector{Vector{Linkinfo{Dim}}}(undef,Dim)

        for μ=1:Dim
            derivatives[μ] = derive_U(line_for_gaugefields,μ)
            numderivatives = length(derivatives[μ] )
            linkinfo_derivatives_right[μ] = Vector{Linkinfo{Dim}}(undef,numderivatives)
            linkinfo_derivatives_leftdag[μ] = Vector{Linkinfo{Dim}}(undef,numderivatives)

            println("$μ ----link")
            show(line_for_gaugefields)
            println("--- derivatives")
            show(derivatives[μ])
            println("\n derivatives")

            for i=1:numderivatives
                rightlinks = get_rightlinks(derivatives[μ][i])
                #shiftedposition = -2*collect(derivatives[μ][i].position)
                #rightlinks = shift_wilsonloop(rightlinks,Tuple(shiftedposition))
                leftlinks = get_leftlinks(derivatives[μ][i])
                #leftlinks = shift_wilsonloop(leftlinks,Tuple(shiftedposition))
                linkinfo_derivatives_right[μ][i] = Linkinfo(rightlinks)
                linkinfo_derivatives_leftdag[μ][i] = Linkinfo(leftlinks')

                println("position")
                println(derivatives[μ][i].position)
                println("")
                println("right")
                show(rightlinks)
                println("")
                println("left")
                show(leftlinks')
                println("")
            end

            println("μ = $μ")

            show(line_for_gaugefields')
            dagline = shift_wilsonloop(line_for_gaugefields',endingpoint)
            show(dagline)
            derivatives_dag[μ] = derive_U(dagline,μ)
            show(derivatives_dag[μ])

            numderivatives = length(derivatives_dag[μ] )
            linkinfo_derivatives_dag_right[μ] = Vector{Linkinfo{Dim}}(undef,numderivatives)
            linkinfo_derivatives_dag_leftdag[μ] = Vector{Linkinfo{Dim}}(undef,numderivatives)
            for i=1:numderivatives
                println("position = ",derivatives_dag[μ][i].position)
                linkinfo_derivatives_dag_right[μ][i] = Linkinfo(get_rightlinks(derivatives_dag[μ][i]))
                linkinfo_derivatives_dag_leftdag[μ][i] = Linkinfo(get_leftlinks(derivatives_dag[μ][i])')
            end
        end

        return new{Dim}( coefficient,matrix_for_spinors,matrix_for_spinors_dagger,
                    matrix_for_spinors_γ5M,
                    line_for_gaugefields,
                    linkinfo,
                    derivatives,
                    linkinfo_derivatives_right,
                    linkinfo_derivatives_leftdag,
                    derivatives_dag,
                    linkinfo_derivatives_dag_right,
                    linkinfo_derivatives_dag_leftdag,
                    startingpoint,
                    endingpoint
                    )
    end

end

function get_positions(d::General_Dirac_term)
    return d.positions
end

function get_directions(d::General_Dirac_term)
    return d.directions
end

function get_isdagvectors(d::General_Dirac_term)
    return d.isdagvectors
end

function get_fermionposition(d::General_Dirac_term)
    return d.endingpoint
end

function get_spinormatrix(d::General_Dirac_term)
    return d.matrix_for_spinors
end

function get_spinormatrix_γ5M(d::General_Dirac_term)
    return d.matrix_for_spinors_γ5M
end

function get_spinormatrix_dagger(d::General_Dirac_term)
    return d.matrix_for_spinors_dagger
end

function get_wilsonline(d::General_Dirac_term)
    return d.line_for_gaugefields
end

function get_coefficient(d::General_Dirac_term)
    return d.coefficient
end

function get_linkinfo(d::General_Dirac_term)
    return d.linkinfo
end

function get_derivatives(d::General_Dirac_term,μ)
    return d.derivatives[μ]
end

function get_linkinfo_derivatives_leftdag(d::General_Dirac_term,μ)
    return d.linkinfo_derivatives_leftdag[μ]
end

function get_linkinfo_derivatives_right(d::General_Dirac_term,μ)
    return d.linkinfo_derivatives_right[μ]
end

function get_derivatives_dag(d::General_Dirac_term,μ)
    return d.derivatives_dag[μ]
end

function get_linkinfo_derivatives_dag_leftdag(d::General_Dirac_term,μ)
    return d.linkinfo_derivatives_dag_leftdag[μ]
end

function get_linkinfo_derivatives_dag_right(d::General_Dirac_term,μ)
    return d.linkinfo_derivatives_dag_right[μ]
end

function Base.show(g::General_Dirac_term{Dim}) where Dim
    println("coefficient: ",g.coefficient)

    println("details: ")
    println("  gauge links: ")
    show(g.line_for_gaugefields)
    println("  matrix for spinors: ")
    display(g.matrix_for_spinors)
    println("\t")
    println("  hermitian-conjugate matrix for spinors: ")
    display(g.matrix_for_spinors_dagger)
    println("\t")
    for μ=1:Dim
        println("dD/dU$μ: ")
        show(g.derivatives[μ])
    end

end

struct Wilson_GeneralDirac_operator{Dim,T,fermion} <: Dirac_operator{Dim}  where T <: AbstractGaugefields
    U::Array{T,1}
    boundarycondition::Vector{Int8}
    _temporary_fermi::Vector{fermion}
    dirac_terms::Vector{General_Dirac_term{Dim}}
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    _temporary_gaugefields::Vector{T}
    verbose_print::Verbose_print
    _temporary_fermion_forCG::Vector{fermion}
end

function Wilson_GeneralDirac_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where {NC,Dim}
    numg = 3
    _temporary_gaugefields = Array{eltype(U),1}(undef,numg)
    for i=1:numg
        _temporary_gaugefields[i] = similar(U[1])
    end

    xtype = typeof(x)
    num = 7
    _temporary_fermi = Array{xtype,1}(undef,num)

    boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])

    inputmode = check_parameters(parameters,"inputmode","easy")

    if inputmode == "expart"
        println("expart mode for constructing Dirac terms")
        dirac_terms = expartmode_constructing_diraclines(parameters,Dim)
    elseif inputmode == "easy"
        println("easy mode for constructing Dirac terms")
        dirac_terms = easymode_constructing_diraclines(parameters,Dim)
    else
        error("inputmode $inputmode is not supported")
    end

    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)
    verbose_level = check_parameters(parameters,"verbose_level",2)
    method_CG = check_parameters(parameters,"method_CG","bicg")

    verbose_print = Verbose_print(verbose_level)

    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    numcg = check_parameters(parameters,"numtempvec_CG",7)
    #numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    return Wilson_GeneralDirac_operator{Dim,eltype(U),xtype}(U,boundarycondition,_temporary_fermi,
        dirac_terms,
        eps_CG,MaxCGstep,
        verbose_level,method_CG,
        _temporary_gaugefields,
        verbose_print,
        _temporary_fermion_forCG
        )
end

function (D::Wilson_GeneralDirac_operator{Dim,T,fermion})(U) where {Dim,T,fermion}
    return Wilson_GeneralDirac_operator{Dim,T,fermion}(U,D.boundarycondition,D._temporary_fermi,
    D.dirac_terms,
    D.eps_CG,
    D.MaxCGstep,
    D.verbose_level,
    D.method_CG,
    D._temporary_gaugefields,
    D.verbose_print,
    D._temporary_fermion_forCG
        )
end

struct DdagD_Wilson_GeneralDirac_operator{Dim,T,fermion} <: DdagD_operator
    dirac::Wilson_GeneralDirac_operator{Dim,T,fermion}
    function DdagD_Wilson_GeneralDirac_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
        return new{Dim,eltype(U),typeof(x)}(Wilson_GeneralDirac_operator(U,x,parameters))
    end

    function DdagD_Wilson_GeneralDirac_operator(D::Wilson_GeneralDirac_operator{Dim,T,fermion}) where {Dim,T,fermion}
        return new{Dim,T,fermion}(D)
    end
end

function get_numterms(w::Wilson_GeneralDirac_operator)
    return length(w.dirac_terms)
end

function Base.getindex(w::Wilson_GeneralDirac_operator,i)
    return w.dirac_terms[i]
end

struct Adjoint_Wilson_GeneralDirac_operator{T} <: Adjoint_Dirac_operator
    parent::T
end

function Base.adjoint(A::T) where T <: Wilson_GeneralDirac_operator
    Adjoint_Wilson_GeneralDirac_operator{typeof(A)}(A)
end

function Base.show(D::Wilson_GeneralDirac_operator{Dim,T,fermion}) where {Dim,T,fermion}
    println("----------------------------------------------")
    println("Wilson Dirac operator: general type")
    println("Method for solving y = D^inv x: ",D.method_CG)
    println("eps for CG: ",D.eps_CG)
    println("Mac steps for CG: ",D.MaxCGstep)
    println("verbose level: ",D.verbose_level)
    numdirac = length(D.dirac_terms)
    println("num. of Dirac terms: ",numdirac)
    for i=1:numdirac
        if i==1
            st ="st"
        elseif i==2
            st ="nd"
        elseif i==3
            st ="rd"
        else
            st ="th"
        end
        println("---------------")
        println("$i-$st term: ")
        show(D.dirac_terms[i])
        println("---------------")
    end
    println("----------------------------------------------")

end

function onelink_diracterm(μ,κ,Dim)
    dirac_terms = General_Dirac_term{Dim}[]

    startingpoint = zeros(Int64,Dim)
    endingpoint = zeros(Int64,Dim)
    endingpoint[μ] = 1

    loop_forward = [(μ,1)]
    line_for_gaugefields = Wilsonline(loop_forward)
    matrix_for_spinors = rminusγ1[:,:,μ]
    dirac_term = General_Dirac_term(κ,matrix_for_spinors,line_for_gaugefields,
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    loop_backward = [(μ,-1)]
    line_for_gaugefields = Wilsonline(loop_backward)
    matrix_for_spinors = rplusγ1[:,:,μ]
    endingpoint[μ] = -1
    dirac_term = General_Dirac_term(κ,matrix_for_spinors,line_for_gaugefields,
                        Tuple(startingpoint),
                        Tuple(endingpoint))
                        push!(dirac_terms,dirac_term)
    return dirac_terms

end

function cloverlink_diracterm_debug(μ,ν,cSW,Dim)
    dirac_terms = General_Dirac_term{Dim}[]
    c = -cSW/(16)

    startingpoint = zeros(Int64,Dim)
    endingpoint = zeros(Int64,Dim)

    γμ = γmatrices[:,:,μ]
    γν = γmatrices[:,:,ν]

    σμν =  (γμ*γν .- γν*γμ )*(1/2)

    loop_righttop = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)])
    loop_lefttop = Wilsonline([(ν,1),(μ,-1),(ν,-1),(μ,1)])
    loop_rightbottom = Wilsonline([(ν,-1),(μ,1),(ν,1),(μ,-1)])
    loop_leftbottom= Wilsonline([(μ,-1),(ν,-1),(μ,1),(ν,1)])

    #loopterm = loop_righttop #correct
    #loopterm = loop_lefttop #not correct
    #loopterm = loop_rightbottom #not correct
    loopterm = loop_leftbottom #not correct
    #loopterm = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)])

    dirac_term = General_Dirac_term(c,σμν,loopterm,
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(-c,σμν,loopterm',
    Tuple(startingpoint),
    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    return dirac_terms

end

function cloverlink_diracterm(μ,ν,cSW,Dim)
    dirac_terms = General_Dirac_term{Dim}[]
    c = -cSW/(16)

    startingpoint = zeros(Int64,Dim)
    endingpoint = zeros(Int64,Dim)

    γμ = γmatrices[:,:,μ]
    γν = γmatrices[:,:,ν]

    σμν =  (γμ*γν .- γν*γμ )*(1/2)

    loop_righttop = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)])
    loop_lefttop = Wilsonline([(ν,1),(μ,-1),(ν,-1),(μ,1)])
    loop_rightbottom = Wilsonline([(ν,-1),(μ,1),(ν,1),(μ,-1)])
    loop_leftbottom= Wilsonline([(μ,-1),(ν,-1),(μ,1),(ν,1)])

    dirac_term = General_Dirac_term(c,σμν,loop_righttop,
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(c,σμν,loop_lefttop,
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(c,σμν,loop_rightbottom,
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(c,σμν,loop_leftbottom,
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(-c,σμν,loop_righttop',
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(-c,σμν,loop_lefttop',
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(-c,σμν,loop_rightbottom',
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    dirac_term = General_Dirac_term(-c,σμν,loop_leftbottom',
                    Tuple(startingpoint),
                    Tuple(endingpoint))

    push!(dirac_terms,dirac_term)

    return dirac_terms

end

function easymode_constructing_diraclines(parameters,Dim)
    nameofterms = check_important_parameters(parameters,"nameofterms","[\"one-link\"]")
    @assert typeof(nameofterms) <: Vector "nameofterms is not a vector. Now $nameofterms"
    println(nameofterms)
    numterms = length(nameofterms)

    dirac_terms =General_Dirac_term{Dim}[]

    for i=1:numterms
        term = nameofterms[i]
        if term == "one-link"
            println("Conventional Wilson terms (1+γ) and (1-γ) are added" )
            κ = check_important_parameters(parameters,"κ")
            for μ=1:Dim
                diracterms = onelink_diracterm(μ,κ,Dim)
                append!(dirac_terms,diracterms)
            end

        elseif term == "clover"
            c = check_important_parameters(parameters,"cSW")
            for μ=1:Dim
                for ν=μ:Dim
                    if ν==μ
                        continue
                    end
                    diracterms = cloverlink_diracterm(μ,ν,c,Dim)
                    append!(dirac_terms,diracterms)
                end
            end
            #println("dd ",length(dirac_terms))
        elseif term == "cloverdebug"
            c = check_important_parameters(parameters,"cSW")
            μ=1
            ν=2
            diracterms = cloverlink_diracterm_debug(μ,ν,c,Dim)
            append!(dirac_terms,diracterms)
            #println("dd ",length(dirac_terms))
        else
            error("term $term is not supported")
        end
    end
    #println("dd ",length(dirac_terms))
    return dirac_terms

end

function expartmode_constructing_diraclines(parameters,Dim)
    dirac_terms = check_important_parameters(parameters,"dirac_terms")
    return dirac_terms
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Wilson_GeneralDirac_operator, T3 <:AbstractFermionfields}
    Wx_general!(y,A.U,x,A)
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_Wilson_GeneralDirac_operator, T3 <:AbstractFermionfields}
    Wdagx_general!(y,A.parent.U,x,A.parent)
end

function Wx_general!(xout::T,U::Array{G,1},x::T,A::Wilson_GeneralDirac_operator)  where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields}
    #temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    tempg = A._temporary_gaugefields
    V = tempg[1]
    tempg1 = tempg[2]
    tempg2 = tempg[3]

    clear_fermion!(xout)
    clear_fermion!(temp2)
    numterms = get_numterms(A)

    for i=1:numterms
        diracterm = A[i]
        position  = get_fermionposition(diracterm)
        xshifted = shift_fermion(x,position)
        wilsonline = get_wilsonline(diracterm)

        linkinfo = get_linkinfo(diracterm)

        gaugefield_fermion_mul!(temp1,linkinfo,U,xshifted,A._temporary_fermi)

        #@time gaugefield_fermion_mul!(temp1,wilsonline,U,xshifted,[temp])

        #println(position)

        #wilsonline = get_wilsonline(diracterm)
        #show(wilsonline)

        #evaluate_gaugelinks!(V,wilsonline,U,[tempg1,tempg2])
        #set_wing_U!(V)

        #Vshift = shift_U(V,position)
        #mul!(temp1,V,xshifted)

        S = get_spinormatrix(diracterm)

        mul!(temp1,S)

        coeff = get_coefficient(diracterm)

        add_fermion!(temp2,coeff,temp1)

        #println(xout[1,1,1,1,1,1])
    end
    #error("d")

    add_fermion!(xout,1,x,-1,temp2)
    set_wing_fermion!(xout,A.boundarycondition)

    #display(xout)
    #    exit()
    return
end

function Wdagx_general!(xout::T,U::Array{G,1},x::T,A::Wilson_GeneralDirac_operator)  where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields}
    #temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    tempg = A._temporary_gaugefields
    V = tempg[1]
    tempg1 = tempg[2]
    tempg2 = tempg[3]

    clear_fermion!(temp2)
    clear_fermion!(xout)
    numterms = get_numterms(A)

    for i=1:numterms
        diracterm = A[i]
        position  = get_fermionposition(diracterm)
        xshifted = shift_fermion(x,position)

        #wilsonline = get_wilsonline(diracterm)
        linkinfo = get_linkinfo(diracterm)
        #positions = get_positions(diracterm)
        #directions = get_directions(diracterm)
        #isdagvectors = get_isdagvectors(diracterm)
        gaugefield_fermion_mul!(temp1,linkinfo,U,xshifted,A._temporary_fermi)

        #clear_U!(V)
        #evaluate_gaugelinks!(V,wilsonline,U,[tempg1,tempg2])
        #set_wing_U!(V)

        #Vshift = shift_U(V,position)
        #mul!(temp1,V,xshifted)

        S = get_spinormatrix_dagger(diracterm)

        mul!(temp1,S)

        coeff = get_coefficient(diracterm)

        add_fermion!(temp2,coeff,temp1)
    end

    add_fermion!(xout,1,x,-1,temp2)
    set_wing_fermion!(xout,A.boundarycondition)
    #set_wing_fermion!(xout,A.boundarycondition)
    #error("rr")
    #display(xout)
    #    exit()
    return
end

const zeroposition = [(0,),(0,0),(0,0,0),(0,0,0,0)]

function gaugefield_fermion_mul!(xout,linkinfo,U::Array{T,1},x,
                temps_fermion::Array{<: AbstractFermionfields{NC,Dim},1}) where
                                                {T<: AbstractGaugefields,Dim,NC}
    positions = linkinfo.positions
    directions = linkinfo.directions
    isdagvectors = linkinfo.isdagvectors
    temp1 = temps_fermion[end]
    temp2 = temps_fermion[end-1]

    #glinks = w
    numlinks = length(positions)
    if numlinks == 0
        substitute_fermion!(xout,x)
        #error("no gaugefield found")
        return
    end

    #println(typeof(position))

    if numlinks == 1
        position = positions[numlinks]
        direction = directions[numlinks]
        isUenddag = isdagvectors[numlinks]

        if position == zeroposition[Dim]
            if isUenddag
                mul!(xout,U[direction]',x)
            else
                mul!(xout,U[direction],x)
            end
        else
            Ushift = shift_U(U[direction],position)
            if isUenddag
                mul!(xout,Ushift',x)
            else
                mul!(xout,Ushift,x)
            end
        end
        return
    else
        position = positions[numlinks]
        direction = directions[numlinks]
        isUenddag = isdagvectors[numlinks]

        if position == zeroposition[Dim]
            if isUenddag
                mul!(temp1,U[direction]',x)
            else
                mul!(temp1,U[direction],x)
            end
        else
            Ushift = shift_U(U[direction],position)
            if isUenddag
                mul!(temp1,Ushift',x)
            else
                mul!(temp1,Ushift,x)
            end
        end
    end

    for j=numlinks-1:-1:2
        position = positions[j]
        direction = directions[j]
        isUdag = isdagvectors[j]

        if position == zeroposition[Dim]
            if isUdag
                mul!(temp2,U[direction]',temp1)
            else
                mul!(temp2,U[direction],temp1)
            end
        else
            Ushift = shift_U(U[direction],position)
            if isUdag
                mul!(temp2,Ushift',temp1)
            else
                mul!(temp2,Ushift,temp1)
            end
        end

        temp2,temp1 = temp1,temp2
    end

    position = positions[1]
    direction = directions[1]
    isU1dag = isdagvectors[1]

    if position == zeroposition[Dim]
        if isU1dag
            mul!(xout,U[direction]',temp1)
        else
            mul!(xout,U[direction],temp1)
        end
    else
        Ushift = shift_U(U[direction],position)
        if isU1dag
            mul!(xout,Ushift',temp1)
        else
            mul!(xout,Ushift,temp1)
        end
    end

end


text:
include("../psmsdst.jl")

"""
    solvegreedy(problem)

Solves the given PSMSDST problem using a purely random approach. While jobs
are available to be programmed, it programs any available job on any machine.

Solutions are generated and compared to the best solution that has been found
so far. If no improvement is made to the solution after a number of
``attempts``, the solver stops and the best solution found is returned.
"""
function solverandom(problem::PSMSDSTProblem,attempts::Int)::Vector{Vector{Int}}

    machines = 1:problem.nummachines

    bestsolution = nothing
    bestcost = typemax(Int64)

    # Number of solutions without improvement
    noimprovement = 0

    # Keep trying solutions until the number of attempts is exceeded
    while noimprovement < attempts
        # Create resulting vector for the solution
        solution = [ Vector{Int}() for machineid in machines ]

        # Create a set for the jobs that have not been programmed yet
        availablejobs = Set{Int}(1:problem.numjobs)

        # Program jobs until a complete solution is reached
        while !isempty(availablejobs)

            # Get a random job and machine
            job = rand(availablejobs)
            machine = rand(machines)

            push!(solution[machine],job)
            delete!(availablejobs,job)
        end

        _, _, cost, _ = evaluatesolution(problem,solution)
        noimprovement += 1

        if cost < bestcost
            bestcost = cost
            bestsolution = solution
            noimprovement = 0
        end
    end

    return bestsolution
end

#
#     getsolver()
#
# Returns a pair containing the name of the solver, and a pre-parameterized
# function that only takes the problem as argument.
#
# Internally it also solves a small problem using the solver, ensuring that
# the code for the solver is compiled ahead of time and there is no reduced
# performance during evaluation.
#
# # Examples
# ```julia-repl
# julia> include("solvers/solver.jl")(param1=42.0,param2=64.0)
# "Solver param1=42 param2=64"=>#1 (generic function with 1 method)
# ```

function getsolver(;attempts=20::Int)::Pair{String,Function}
    solver = problem -> solverandom(problem,attempts)

    presolve(solver)

    return "Random n=$attempts"=>solver
end

max_stars_repo_name: mkualquiera/huepl

max_issues_repo_name: mkualquiera/huepl

max_forks_repo_name: mkualquiera/huepl


text:
using Weave

cd(@__DIR__)
weave("backscatter_energy_depth.ipynb", fig_ext=".svg")
weave("backscatter_mc_browning1991.ipynb", fig_ext=".svg")
weave("backscatter_mc_browning1994.ipynb", fig_ext=".svg")
weave("backscatter_mc.ipynb", fig_ext=".svg")
weave("phirhoz_mc.ipynb", fig_ext=".svg")
#weave("visualize_mc.ipynb", fig_ext=".svg")


text:
using EqualityConstrainedMinimization
using LinearAlgebra

#Generates ECM problem with feasible starting point
function Generate_ECM(_n, _p, A_bound)
    _x = rand(_n)
    _A = rand(_p, _n) * A_bound

    it = 0
    #Ensure A has full rank
    while rank(_A) != max(_n, _p) && it < 100 #Cap generation at 100 tries. The odds of needing this many are very low
        _A = rand(_p, _n) * A_bound
        it += 1
    end

    _b = _A * _x

    return _A, _x, _b
end

#Test objective for constrained problem (see Boyd and Vandenberghe problem 10.15)
function Test_Objective_Con(_x)
    _n = length(_x)
    sum = 0

    for i = 1:_n
        sum += _x[i] * log(_x[i])
    end

    return sum
end

l = 5
w = 3 #Constraint coefficient mtx will be an n x m matrix
ECM = Generate_ECM(l, w, 1)
mtx = ECM[1] #the A in Ax = b
x = ones(l)*0.5 #ECM[2] #set equal to ECM[2] to test feasible start
vect = ECM[3] #the b in Ax = b

f(_x) = Test_Objective_Con(_x)

#Solves the problem "minimize f(x) subject to Ax = b." Since b is given, we assume infeasible start
IF_Sol = Solve_ECM(f, x, A = mtx, b = vect, verbose = false)

#Solves the problem "minimize f(x) subject to Ax = b." Since b not given, we assume feasible start
x = ECM[2]
F_Sol = Solve_ECM(f, x, A = mtx, verbose = false)

#Solves the problem "minimize f(x)."
#Setting force_gradient = true will force the gradient method
#Do this if ϵ is relatively large.
#Otherwise, the solver will first use Gradient Descent and then switch to Newton's Method
U_Sol = Solve_ECM(f, x, verbose = false)

max_stars_repo_name: ihawn/EqualityConstrainedMinimization

max_issues_repo_name: ihawn/EqualityConstrainedMinimization

max_forks_repo_name: ihawn/EqualityConstrainedMinimization


text:
# usage from repository root:
#
# julia --project=. examples/fasync_tumour_invasion.jl
#

using BooleanNetworks

import Random
Random.seed!(1234)

model = "models/Tumour_invasion.bnet"
outputs = ["Apoptosis"; "CellCycleArrest"; "Invasion"; "Metastasis"; "Migration"]
free_nodes = ["DNAdamage", "ECMicroenv"]

bn = load_bnet("models/Tumour_invasion.bnet")

outputs = [bn.index[node] for node in outputs]
free_nodes = [bn.index[node] for node in free_nodes]

x0 = zerocfg(bn)
for i in free_nodes
    x0[i] = true
end

# warm-up
fasync_simulations(bn, outputs, 1, 10, x0)

#using BenchmarkTools

nb_sims = 100000
maxsteps = 500
for _ in 1:4
    println("Performing $nb_sims simulations of at most $maxsteps from fixed initial configuration")
    @time fasync_simulations(bn, outputs, nb_sims, maxsteps, x0);
end

using Profile
@profile fasync_simulations(bn, outputs, nb_sims, maxsteps, x0);
Profile.print()

result = fasync_simulations(bn, outputs, nb_sims, maxsteps, x0);
ratios = outputs_ratios(result, outputs, bn)
println(ratios)


text:
# Constructor for MLEControl struct

"""
    MLESettings(y, model, init, optimizer, ci)
Wrapper function to specify estimation settings.

* `y`: Time series
* `model`: Model specification
* `init`: Initial values (vector or parameter).
* `optimizer`: Optimization routine. "BFGS", "LBFGS" or "NelderMead"
* `ci`: Indicator: Shall confidence intervals be computed?
* `maxEval`: Maximum number of likelihood evaluations

# Example
```julia-repl
MLESettings(y, model, ci = true, maxEval = 1e9)
```

If the argument `init` is not given, valid initial values are chosen.
See also [MLEControl](@ref).
"""
function MLESettings(y::Array{T, 1} where T<:Integer,
    model::INGARCHModel,
    init::Array{T, 1} where T<:AbstractFloat = Array{Float64, 1}([]);
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    p = length(model.pastObs)
    q = length(model.pastMean)
    r = length(model.external)
    nb = model.distr == "NegativeBinomial"

    nPar = 1 + p + q + r + nb + model.zi

    if (length(init) != nPar) & (length(init) > 0)
        println("Number of initial values does not match number of parameters.")
        println("Switched to default initial values.")
        init = Array{Float64, 1}(undef, 0)
    end

    if length(init) == 0
        β0init = mean(y)/2
        if model.link == "Log"
            β0init = log(β0init)
        end

        M = p + q
        if M > 0
            temp = [collect(p:-1:1); collect(q:-1:1)]
            temp = temp./(2*sum(temp))
        else
            temp = Array{Float64, 1}([])
        end

        init = [β0init; temp; fill(0.05, r); fill(5.0, nb); ifelse(model.zi, 0.1, Array{Float64, 1}([]))]
    end

    init = θ2par(init, model)

    if !(optimizer in ["BFGS", "LBFGS", "NelderMead"])
        println("Invalid optimizing method. Set to BFGS.")
        optimizer = "BFGS"
    end

    MLEControl(init, optimizer, ci, Int64(round(abs(maxEval))))
end

function MLESettings(y::Array{T, 1} where T<:Integer,
    model::INARCHModel,
    init::Array{T, 1} where T<:AbstractFloat = Array{Float64, 1}([]);
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    p = length(model.pastObs)
    r = length(model.external)
    nb = model.distr == "NegativeBinomial"

    nPar = 1 + p + r + nb + model.zi

    if (length(init) != nPar) & (length(init) > 0)
        println("Number of initial values does not match number of parameters.")
        println("Switched to default initial values.")
        init = Array{Float64, 1}(undef, 0)
    end

    if length(init) == 0
        β0init = mean(y)/2
        if model.link == "Log"
            β0init = log(β0init)
        end

        if p > 0
            temp = collect(p:-1:1)
            temp = temp./(2*sum(temp))
        else
            temp = Array{Float64, 1}([])
        end

        init = [β0init; temp; fill(0.05, r); fill(5.0, nb); ifelse(model.zi, 0.1, Array{Float64, 1}([]))]
    end

    init = θ2par(init, model)

    if !(optimizer in ["BFGS", "LBFGS", "NelderMead"])
        println("Invalid optimizing method. Set to BFGS.")
        optimizer = "BFGS"
    end

    MLEControl(init, optimizer, ci, Int64(round(abs(maxEval))))
end

function MLESettings(y::Array{T, 1} where T<:Integer,
    model::IIDModel,
    init::Array{T, 1} where T<:AbstractFloat = Array{Float64, 1}([]);
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    r = length(model.external)
    nb = model.distr == "NegativeBinomial"

    nPar = 1 + r + nb + model.zi

    if (length(init) != nPar) & (length(init) > 0)
        println("Number of initial values does not match number of parameters.")
        println("Switched to default initial values.")
        init = Array{Float64, 1}(undef, 0)
    end

    if length(init) == 0
        β0init = mean(y)/2
        if model.link == "Log"
            β0init = log(β0init)
        end

        init = [β0init; fill(0.05, r); fill(5.0, nb); ifelse(model.zi, 0.1, Array{Float64, 1}([]))]
    end

    init = θ2par(init, model)

    if !(optimizer in ["BFGS", "LBFGS", "NelderMead"])
        println("Invalid optimizing method. Set to BFGS.")
        optimizer = "BFGS"
    end

    MLEControl(init, optimizer, ci, Int64(round(abs(maxEval))))
end

function MLESettings(y::Array{T, 1} where T<:Integer,
    model::INARMAModel,
    init::Array{T, 1} where T<:AbstractFloat = Array{Float64, 1}([]);
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    p = length(model.pastObs)
    q = length(model.pastMean)
    r = length(model.external)
    nϕ = sum(model.distr .== "NegativeBinomial")
    if nϕ == 2
        if r == 0
            nϕ -= 1
        else
            if sum(model.external) == 0
                nϕ -= 1
            end
        end
    end

    nPar = 1 + p + q + r + nϕ + model.zi

    if (length(init) != nPar) & (length(init) > 0)
        println("Number of initial values does not match number of parameters.")
        println("Switched to default initial values.")
        init = Array{Float64, 1}(undef, 0)
    end

    if length(init) == 0
        β0init = mean(y)/2
        if model.link[1] == "Log"
            β0init = log(β0init)
        end

        M = p + q
        if M > 0
            temp = [collect(p:-1:1); collect(q:-1:1)]
            temp = temp./(2*sum(temp))
        else
            temp = Array{Float64, 1}([])
        end

        init = [β0init; temp; fill(0.05, r); fill(5.0, nϕ); ifelse(model.zi, 0.1, Array{Float64, 1}([]))]
    end

    init = θ2par(init, model)

    if !(optimizer in ["BFGS", "LBFGS", "NelderMead"])
        println("Invalid optimizing method. Set to BFGS.")
        optimizer = "BFGS"
    end

    MLEControl(init, optimizer, ci, Int64(round(abs(maxEval))))
end

function MLESettings(y::Array{T, 1} where T<:Integer,
    model::INARModel,
    init::Array{T, 1} where T<:AbstractFloat = Array{Float64, 1}([]);
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    p = length(model.pastObs)
    r = length(model.external)
    nϕ = sum(model.distr .== "NegativeBinomial")
    if nϕ == 2
        if r == 0
            nϕ -= 1
        else
            if sum(model.external) == 0
                nϕ -= 1
            end
        end
    end

    nPar = 1 + p + r + nϕ + model.zi

    if (length(init) != nPar) & (length(init) > 0)
        println("Number of initial values does not match number of parameters.")
        println("Switched to default initial values.")
        init = Array{Float64, 1}(undef, 0)
    end

    if length(init) == 0
        β0init = mean(y)/2
        if model.link[1] == "Log"
            β0init = log(β0init)
        end

        if p > 0
            temp = collect(p:-1:1)
            temp = temp./(2*sum(temp))
        else
            temp = Array{Float64, 1}([])
        end

        init = [β0init; temp; fill(0.05, r); fill(5.0, nϕ); ifelse(model.zi, 0.1, Array{Float64, 1}([]))]
    end

    init = θ2par(init, model)

    if !(optimizer in ["BFGS", "LBFGS", "NelderMead"])
        println("Invalid optimizing method. Set to BFGS.")
        optimizer = "BFGS"
    end

    MLEControl(init, optimizer, ci, Int64(round(abs(maxEval))))
end

function MLESettings(y::Array{T, 1} where T<:Integer,
    model::INMAModel,
    init::Array{T, 1} where T<:AbstractFloat = Array{Float64, 1}([]);
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    q = length(model.pastMean)
    r = length(model.external)
    nϕ = sum(model.distr .== "NegativeBinomial")
    if nϕ == 2
        if r == 0
            nϕ -= 1
        else
            if sum(model.external) == 0
                nϕ -= 1
            end
        end
    end

    nPar = 1 + q + r + nϕ + model.zi

    if (length(init) != nPar) & (length(init) > 0)
        println("Number of initial values does not match number of parameters.")
        println("Switched to default initial values.")
        init = Array{Float64, 1}(undef, 0)
    end

    if length(init) == 0
        β0init = mean(y)/2
        if model.link[1] == "Log"
            β0init = log(β0init)
        end

        if q > 0
            temp = collect(q:-1:1)
            temp = temp./(2*sum(temp))
        else
            temp = Array{Float64, 1}([])
        end

        init = [β0init; temp; fill(0.05, r); fill(5.0, nϕ); ifelse(model.zi, 0.1, Array{Float64, 1}([]))]
    end

    init = θ2par(init, model)

    if !(optimizer in ["BFGS", "LBFGS", "NelderMead"])
        println("Invalid optimizing method. Set to BFGS.")
        optimizer = "BFGS"
    end

    MLEControl(init, optimizer, ci, Int64(round(abs(maxEval))))
end

function MLESettings(y::Array{T, 1} where T<:Integer,
    model::T where T<:CountModel,
    init::parameter;
    optimizer::String = "BFGS",
    ci::Bool = false,
    maxEval::T where T<:AbstractFloat = 1e10)

    initNew = par2θ(init, model)
    MLESettings(y, model, initNew, optimizer = optimizer, ci = ci, maxEval = maxEval)
end


text:
"""
    without_incoming(G::AbstractGraph, X::Set)

Return a graph where all the arrows pointing to nodes in the set `X` in graph `G` have been
removed.
"""
function without_incoming(G::AbstractGraph, X::Set)
    E = collect(LightGraphs.edges(G))
    filter!(e -> !(e.dst in X), E)
    return SimpleDiGraph(E)
end

"""
    without_outgoing(G::AbstractGraph, X::Set)

Return a graph where all the arrows coming from nodes in the set `X` in graph `G` have been
removed.
"""
function without_outgoing(G::AbstractGraph, X::Set)
    E = collect(LightGraphs.edges(G))
    filter!(e -> !(e.src in X), E)
    return SimpleDiGraph(E)
end

"""
    d_separated(G::AbstractGraph, X::Int, Y::Int, Z::Set)::Bool

Return whether the elements in `Z` block all paths from `X` to `Y` in graph `G`.
In other words, whether `X` and `Y` are _d_-separated by `Z`, normally written as
`X ⊥⊥ Y ¦ Z`.

Specifically, return whether for all bi-directional paths `path` from `X` to `Y`
[^Bareinboim2015]:

1. `path` contains at least one arrow-emitting node that is in `Z` _or_
2. `path` contains at least one collision node that is outside `Z` and has no descendant
    in `Z`.

[^Bareinboim2015]: Bareinboim, E., & Pearl, J. (2016). Causal inference and the data-fusion problem. Proceedings of the National Academy of Sciences, 113(27), 7345-7352.

"""
function d_separated(G::AbstractGraph, X::Int, Y::Int, Z::Set; verbose=false)::Bool
    @assert isdisjoint(X, Y)
    @assert isdisjoint(X, Z)
    @assert isdisjoint(Y, Z)
    @assert length(nodes(G)) != 0
    return CausalInference.dsep(G, X, Y, Z; verbose)
end

function d_separated(G::AbstractGraph, X::Int, Y::Int, Z::Vector; verbose=false)::Bool
    return d_separated(G, X, Y, Set(Z); verbose)
end

function d_separated(G::AbstractGraph, X, Y, Z::Set; verbose=false)
    if length(X) == 1
        X = first(X)
    else
        error("Not implemented")
    end
    if length(Y) == 1
        Y = first(Y)
    else
        error("Not implemented")
    end
    d_separated(G, X, Y, Z; verbose)
end

function nodes(G::AbstractGraph)
    E = collect(edges(G))
    nodes = unique(union(getproperty.(E, :src), getproperty.(E, :dst)))
    return nodes
end

"""
    graph(edges::Vector{Pair{Symbol, Symbol}})::MetaDiGraph

Return a MetaDiGraph from the `edges` where all the vertices have a `:name` property.
See also `name2node`[@ref].

# Example
```jldoctest
julia> edges = [:U => :V, :V => :W];

julia> G = graph(edges);

julia> get_prop(G, 1, :name)
:U

julia> name2node(G, :U)
1
"""
function graph(edges::Vector{Pair{Symbol, Symbol}})::MetaDiGraph
    N = unique(union(first.(edges), last.(edges)))
    sym2int = Dict(zip(N, 1:length(N)))
    edges = [Edge(sym2int[src] => sym2int[dst]) for (src, dst) in edges]
    G = SimpleDiGraph(edges)
    G = MetaDiGraph(G)

    int2sym = Dict(zip(1:length(N), N))
    foreach(k -> set_prop!(G, k, :name, int2sym[k]), keys(int2sym))

    return G
end

"""
    name2node(G::MetaDiGraph, name::Symbol)::Int

Return the node for a given `name`.
"""
function name2node(G::MetaDiGraph, name::Symbol)::Int
    V = filter_vertices(G, :name, name)
    V = collect(V)
    @assert length(V) == 1
    return first(V)
end


text:
####################################################################################################
####################################################################################################
#Simple Pointwise Operators
exp(dom::AbstractDomain) = PointwiseOperator(dom, [exp], [exp])
expm1(dom::AbstractDomain) = PointwiseOperator(dom, [expm1], [exp])
log(dom) = PointwiseOperator(dom, [log], [inv])
log1p(dom) = PointwiseOperator(dom, [log1p], [x -> 1/(1 + x)])
sin(dom) = PointwiseOperator(dom, [cos], [sin])
cos(dom) = PointwiseOperator(dom, [cos], [sin, -])

add(dom, val::Union{Number, AbstractArray{Number}}) = PointwiseOperator(dom, [x -> +(x, val)], [identity])
scale(dom, val::Union{Number, AbstractArray{Number}}) = PointwiseOperator(dom, [x -> *(x, val)], [x -> *(x, val)])

####################################################################################################
####################################################################################################
#Simple Linear Operators
struct HartleyOperator <: LinearOperator
    domain::StructuredDomain{N, true} where N
end
target(op::HartleyOperator) = getcodomain(op.domain)
apply(op::HartleyOperator, val) = hartley(val)
adjoint(op::HartleyOperator, val) = apply(op, val)

hartley(op::AbstractOperator) = combine_operators(HartleyOperator(target(op)), op)

struct FourierOperator <: LinearOperator
    domain::StructuredDomain{N, true} where N
end
target(op::FourierOperator) = getcodomain(op.domain)
apply(op::FourierOperator, val) = fft(val)
#adjoint(op::FourierOperator, val) = fft(conj.(val))

####################################################################################################
####################################################################################################
#Multi-Field Operators
struct AdditionOperator <: AbstractOperator
    domain
    target
    operators
end

struct SubtractionOperator <: AbstractOperator
    domain
    target
    operators
end

struct MultiplicationOperator <: AbstractOperator
    domain
    target
    operators
end

struct DivisionOperator <: AbstractOperator
    domain
    target
    operators
end

function +(op1::AbstractOperator, op2::AbstractOperator)
    @assert target(op1) == target(op2)
    dom = (domain(op1), domain(op2))
    return AdditionOperator(dom, target(op1), (op1, op2))
end

function -(op1::AbstractOperator, op2::AbstractOperator)
    @assert target(op1) == target(op2)
    dom = (domain(op1), domain(op2))
    return SubtractionOperator(dom, target(op1), (op1, op2))
end

function *(op1::AbstractOperator, op2::AbstractOperator)
    @assert target(op1) == target(op2)
    dom = (domain(op1), domain(op2))
    return MultiplicationOperator(dom, target(op1), (op1, op2))
end

function /(op1::AbstractOperator, op2::AbstractOperator)
    @assert target(op1) == target(op2)
    dom = (domain(op1), domain(op2))
    return DivisionOperator(dom, target(op1), (op1, op2))
end

function apply(op::AdditionOperator, vals)
    #TODO replace zeros(size(target)) with zeros(target)
    res = zeros(size(target(op)))
    for (operation, val) in zip(op.operators, vals)
        res .+= apply(operation, val)
    end
    return res
end

function apply(op::SubtractionOperator, (val1, val2))
    return apply(op.operators[1], val1) .- apply(op.operators[2], val2)
end

function apply(op::MultiplicationOperator, vals)
    #TODO replace zeros(size(target)) with zeros(target)
    res = ones(size(target(op)))
    for (operation, val) in zip(op.operators, vals)
        res .*= apply(operation, val)
    end
    return res
end

function apply(op::SubtractionOperator, (val1, val2))
    return apply(op.operators[1], val1) ./ apply(op.operators[2], val2)
end
####################################################################################################
####################################################################################################
#Convenience Methods
for f in [:exp, :expm1, :log, :log1p, :sin, :cos]
    @eval $f(op::AbstractOperator) = combine_operators($f(target(op)), op)
end


text:
module TimeZones

using Dates
using Printf
using Serialization
using RecipesBase: RecipesBase, @recipe
using Unicode

import Dates: TimeZone, UTC

export TimeZone, @tz_str, istimezone, FixedTimeZone, VariableTimeZone, ZonedDateTime,
    DateTime, Date, Time, UTC, Local, TimeError, AmbiguousTimeError, NonExistentTimeError,
    UnhandledTimeError,
    # discovery.jl
    timezone_names, all_timezones, timezones_from_abbr, timezone_abbrs,
    next_transition_instant, show_next_transition,
    # accessors.jl
    timezone, hour, minute, second, millisecond,
    # adjusters.jl
    firstdayofweek, lastdayofweek,
    firstdayofmonth, lastdayofmonth,
    firstdayofyear, lastdayofyear,
    firstdayofquarter, lastdayofquarter,
    # Re-export from Dates
    yearmonthday, yearmonth, monthday, year, month, week, day, dayofmonth,
    # conversion.jl
    now, today, todayat, astimezone,
    # local.jl
    localzone,
    # ranges.jl
    guess

const PKG_DIR = dirname(@__DIR__)
const DEPS_DIR = joinpath(PKG_DIR, "deps")

# TimeZone types used to disambiguate the context of a DateTime
# abstract type UTC <: TimeZone end  # Already defined in the Dates stdlib
abstract type Local <: TimeZone end

function __init__()
    # Initialize the thread-local TimeZone cache (issue #342)
    _reset_tz_cache()

    # Base extension needs to happen everytime the module is loaded (issue #24)
    Dates.CONVERSION_SPECIFIERS['z'] = TimeZone
    Dates.CONVERSION_SPECIFIERS['Z'] = TimeZone
    Dates.CONVERSION_DEFAULTS[TimeZone] = ""
    Dates.CONVERSION_TRANSLATIONS[ZonedDateTime] = (
        Year, Month, Day, Hour, Minute, Second, Millisecond, TimeZone,
    )

    global ISOZonedDateTimeFormat = DateFormat("yyyy-mm-ddTHH:MM:SS.ssszzz")
end

include("compat.jl")
include("utils.jl")
include("indexable_generator.jl")

include("class.jl")
include("utcoffset.jl")
include(joinpath("types", "timezone.jl"))
include(joinpath("types", "fixedtimezone.jl"))
include(joinpath("types", "variabletimezone.jl"))
include(joinpath("types", "zoneddatetime.jl"))
include("exceptions.jl")
include(joinpath("tzdata", "TZData.jl"))
Sys.iswindows() && include(joinpath("winzone", "WindowsTimeZoneIDs.jl"))
include("build.jl")
include("interpret.jl")
include("accessors.jl")
include("arithmetic.jl")
include("io.jl")
include("tzfile.jl")
include("adjusters.jl")
include("conversions.jl")
include("local.jl")
include("ranges.jl")
include("discovery.jl")
include("rounding.jl")
include("parse.jl")
include("plotting.jl")
include("deprecated.jl")

end # module


text:
using FractalTools
using Plots

# Define function
f(x) = 4 - (x - 2)^2 + 1

# Interpolation points
xi, dx, xf = 0, 1, 4
x = collect(xi : dx : xf)
y = f.(x)
d = 0.5 * ones(length(x) - 1)

# Construct interpolation
interp = interpolate(x, y, d)
ws = interp.ifs.ws

# End points
xc = [0, 4, 4, 0, 0]
yc = [0, 0, 4, 4, 0]
pc = [xc yc]

for niter in 1 : niters
    xt = collect(xi : 0.1dx : xf)
    plt = plot(xt, f.(xt))
    pc = vcat([mapslices(wt, pc, dims=2) for wt in ws]...)
end

function transform_plot()
    # plt = plot(getindex.(pc, 1), getindex.(pc,2), marker=(:circle, 3))
    # xt = collect(xi : 0.1dx : xf)
    # plot!(xt, f.(xt))
    for wt in ws
        wp = []
        for (xi, yi) in zip(xc, yc)
            push!(wp, wt([xi, yi]))
        end
        plot!(getindex.(wp, 1), getindex.(wp,2), marker=(:circle, 3), subplot=1)
    end
    plt
end

transform_plot()


text:
module Layouts

    import ..Log4jl: Layout, header, footer,
                     StringLayout,
                     Event, message, level, timestamp, fqmn, logger, marker,
                     Message, formatted,
                     LOG4JL_LINE_SEPARATOR, getbacktrace

    import Base: serialize, string, show

    include("layouts/basic.jl")
    include("layouts/pattern.jl")
    include("layouts/serialized.jl")

end


text:
#=
decompose a matrix C into a sum of a sparse binary matrix A and a matrix B with small
spectral norm, constraining the number of nonzeros in A

given matrix C, solve:
min specnorm(C - A) :
Aᵢⱼ ∈ {0,1}
sum(A) = k
=#

struct MatrixDecomposition <: ExampleInstance
    m::Int
    n::Int
    use_nat::Bool # use nuclear norm cone, else PSD cone formulation
end

function build(inst::MatrixDecomposition)
    (m, n) = (inst.m, inst.n)
    B_rank = 1 + div(m, 2)
    @assert 1 <= B_rank <= m <= n

    # generate data
    A_sparsity = min(1.0, 2 / n)
    A0 = sprand(Bool, m, n, A_sparsity)
    k = nnz(A0)
    B0 = randn(m, B_rank) * randn(B_rank, n)
    C = A0 + B0

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, A[1:m, 1:n], Bin)
    JuMP.@constraint(model, sum(A) == k)

    JuMP.@variable(model, u)
    JuMP.@objective(model, Min, u)

    add_nonsymm_specnuc(inst.use_nat, false, u, C - A, model)

    # save for use in tests
    model.ext[:u] = u
    model.ext[:A] = A
    model.ext[:B] = C - A

    return model
end

function test_extra(inst::MatrixDecomposition, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    A = JuMP.value.(model.ext[:A])
    @test A ≈ round.(Int, A) atol = tol rtol = tol
    @test all(-tol .<= A .<= 1 + tol)
    # check nuclear norm value
    tol = eps()^0.2
    u = JuMP.value.(model.ext[:u])
    B = JuMP.value.(model.ext[:B])
    @test u ≈ maximum(svdvals(B)) atol = tol rtol = tol
    @test JuMP.objective_value(model) ≈ u
    return
end


text:
function get_graph(graphname, size; degree, filling, seed)
    if graphname=="diag"
        return DiagGraphConfig(; filling, n=size, m=size, seed=seed)
    elseif graphname == "regular"
        return RegularGraphConfig(; degree, size=size, seed=seed)
    elseif graphname == "square"
        return SquareGraphConfig(; filling, m=size, n=size, seed=seed)
    elseif graphname == "misproject"
        return MISProjectGraphConfig(; n=size, index=seed)
    elseif graphname == "mapped-regular"
        return MappedRegularGraphConfig(; degree, size, seed)
    else
        return SmallGraphConfig(graphname)
    end
end

function unified_solve(nv::Int, edges::Vector{Vector{Int}}, problem::String, property::String;
            weights, openvertices, fixedvertices,
            # contraction order optimizer
            sc_target, sc_weight, ntrials, niters, rw_weight, nslices, βs=0.01:0.05:30,
            cudadevice::Int,
        )
    graph = SimpleGraph(nv)
    for (i, j) in edges
        add_edge!(graph, i, j)
    end
    PT = parseproblem(problem)
    weights = weights === nothing ? NoWeight() : weights
    instance = PT(graph; weights, openvertices=openvertices, fixedvertices=fixedvertices,
        optimizer=TreeSA(; sc_target, sc_weight, ntrials, niters, rw_weight, nslices, βs),
        simplifier=MergeGreedy()
    )

    property = GraphUtilities.parseproperty(property)

    println("time, space, RW complexitys are $(timespacereadwrite_complexity(instance))")
    if cudadevice >=0
        CUDA.device!(cudadevice)
    end
    return solve(instance, property; usecuda=cudadevice>=0)[]
end

function loadjson_graph(dict)
    nv = dict["nv"]
    edges = dict["edges"]
    graph = SimpleGraph(nv)
    for (i, j) in edges
        add_edge!(graph, i, j)
    end
    return graph
end

function loadjson_optimizer(dict)
    method = dict["method"]
    dict = dict["method.$method"]  # extract the related method
    if method == "TreeSA"
        return TreeSA(;
            sc_target=dict["sc_target"],
            sc_weight=dict["sc_weight"],
            rw_weight=dict["rw_weight"],
            ntrials=dict["ntrials"],
            niters=dict["niters"],
            nslices=dict["nslices"],
            βs=dict["betas"]
        )
    elseif method == "GreedyMethod"
        return GreedyMethod(;
            nrepeat=dict["nrepeat"]
        )
    else
        error("")
    end
end

function tojson(data)
    if data isa Tropical
        return Dict("size"=>data.n)
    elseif data isa ExtendedTropical{K,<:Tropical} where K
        return [Dict("size"=>x.n) for x in data.orders]
    elseif data isa Real
        return Dict("count"=>data)
    elseif data isa CountingTropical{T, <:Real} where T
        return Dict("size"=>data.n, "count"=>data.c)
    elseif data isa TruncatedPoly{K, <:Real} where {K}
        K = length(data.coeffs)
        return [Dict("size"=>size, "count"=>count) for (size, count) in zip(data.maxorder-K+1:data.maxorder, data.coeffs)]
    elseif data isa Polynomial
        return [Dict("size"=>i-1, "count"=>c) for (i, c) in enumerate(data.coeffs)]
    elseif data isa CountingTropical{T, <:ConfigSampler} where T
        return Dict("size"=>data.n, "config"=>[data.c.data...])
    elseif data isa ExtendedTropical{K, <:CountingTropical{T, <:ConfigSampler}} where {K,T}
        return [Dict(
                    "size"=>di.n,
                    "config"=>[di.c.data...]
                    ) for di in data.orders
                ]
    elseif data isa CountingTropical{T, <:ConfigEnumerator} where T
        return Dict("size"=>data.n, "configs"=>[[x...] for x in data.c])
    elseif data isa TruncatedPoly{K, <:ConfigEnumerator} where K
        sizes = collect(data.maxorder-length(data.coeffs)+1:data.maxorder)
        configs = [data.coeffs...]
        return [Dict("size"=>size, "configs"=>[[x...] for x in config]) for (size, config) in zip(sizes, configs)]
    elseif data isa ConfigEnumerator
        return Dict("configs"=>[[x...] for x in data])
    elseif data isa CountingTropical{T, <:SumProductTree} where T
        return Dict("size"=>data.n, "configs"=>[[x...] for x in collect(data.c)])
    elseif data isa TruncatedPoly{K, <:SumProductTree} where K
        sizes = collect(data.maxorder-length(data.coeffs)+1:data.maxorder)
        configs = collect.([data.coeffs...])
        return [Dict("size"=>size, "configs"=>[[x...] for x in config]) for (size, config) in zip(sizes, configs)]
    elseif data isa SumProductTree
        return Dict("configs"=>[[x...] for x in collect(data)])
    else
        error("unknown data type: `$(typeof(data))`.")
    end
end

function application(dict)
    dict = verify(application_verifier, dict)
    api = dict["api"]
    if api == "solve"
        return app_solve(dict["api.solve"])
    elseif api == "graph"
        return app_graph(dict["api.graph"])
    elseif api == "opteinsum"
        return app_opteinsum(dict["api.opteinsum"])
    elseif api == "help"
        Dict("help"=>string(application_verifier))
    end
end

function app_opteinsum(dict)
    code = GenericTensorNetworks.OMEinsum.DynamicEinCode(dict["inputs"], dict["output"])
    optimizer = loadjson_optimizer(dict["optimizer"])
    simp = get(dict, "simplifier", "")  # optional
    simplifier = if simp == "MergeGreedy"
        MergeGreedy()
    elseif simp == "MergeVectors"
        MergeVectors()
    else
        nothing
    end
    GenericTensorNetworks.OMEinsumContractionOrders._todict(optimize_code(code, dict["sizes"], optimizer, simplifier))
end

function app_solve(dict)
    graph = loadjson_graph(dict["graph"])
    optimizer = loadjson_optimizer(dict["optimizer"])
    PROB = parseproblem(dict["problem"])
    property = parseproperty(dict["property"])

    instance = PROB(graph;
            weights=get(dict, "weights", NoWeight()),  # weights is optional
            openvertices=dict["openvertices"],
            fixedvertices=dict["fixedvertices"],
            optimizer,
            simplifier=MergeGreedy()
        )
    println("time, space, RW complexitys are $(timespacereadwrite_complexity(instance))")

    cudadevice = dict["cudadevice"]
    cudadevice >= 0 && CUDA.device!(cudadevice)
    return tojson(solve(instance, property; usecuda=cudadevice>=0)[])
end

function app_graph(dict)
    type = dict["type"]
    cfg = dict["type.$type"]
    g = if type == "kings"
        Random.seed!(cfg["seed"])
        random_diagonal_coupled_graph(cfg["m"], cfg["n"], cfg["filling"])
    elseif type == "square"
        Random.seed!(cfg["seed"])
        random_square_lattice_graph(cfg["m"], cfg["n"], cfg["filling"])
    elseif type == "smallgraph"
        Random.seed!(cfg["seed"])
        smallgraph(Symbol(cfg["name"]))
    elseif type == "regular"
        random_regular_graph(cfg["n"], cfg["d"], seed=cfg["seed"])
    else
        error("unkown graph type: $type")
    end
    return Dict("nv"=>nv(g), "edges"=>[[e.src, e.dst] for e in edges(g)])
end


text:
export iauDtf2d
"""
Encode date and time fields into 2-part Julian Date (or in the case
of UTC a quasi-JD form that includes special provision for leap
seconds).

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   scale     char[]  time scale ID (Note 1)
   iy,im,id  int     year, month, day in Gregorian calendar (Note 2)
   ihr,imn   int     hour, minute
   sec       double  seconds

Returned:
   d1,d2     double  2-part Julian Date (Notes 3,4)

Returned (function value):
            int     status: +3 = both of next two
                              +2 = time is after end of day (Note 5)
                              +1 = dubious year (Note 6)
                              0 = OK
                              -1 = bad year
                              -2 = bad month
                              -3 = bad day
                              -4 = bad hour
                              -5 = bad minute
                              -6 = bad second (<0)

Notes:

   1. scale identifies the time scale.  Only the value "UTC" (in upper
   case) is significant, and enables handling of leap seconds (see
   Note 4).

   2. For calendar conventions and limitations, see iauCal2jd.

   3. The sum of the results, d1+d2, is Julian Date, where normally d1
   is the Julian Day Number and d2 is the fraction of a day.  In the
   case of UTC, where the use of JD is problematical, special
   conventions apply:  see the next note.

   4. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The SOFA internal convention is that
   the quasi-JD day represents UTC days whether the length is 86399,
   86400 or 86401 SI seconds.  In the 1960-1972 era there were
   smaller jumps (in either direction) each time the linear UTC(TAI)
   expression was changed, and these "mini-leaps" are also included
   in the SOFA convention.

   5. The warning status "time is after end of day" usually means that
   the sec argument is greater than 60.0.  However, in a day ending
   in a leap second the limit changes to 61.0 (or 59.0 in the case
   of a negative leap second).

   6. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See iauDat for further details.

   7. Only in the case of continuous and regular time scales (TAI, TT,
   TCG, TCB and TDB) is the result d1+d2 a Julian Date, strictly
   speaking.  In the other cases (UT1 and UTC) the result must be
   used with circumspection;  in particular the difference between
   two such results cannot be interpreted as a precise time
   interval.

Called:
   iauCal2jd    Gregorian calendar to JD
   iauDat       delta(AT) = TAI-UTC
   iauJd2cal    JD to Gregorian calendar

This revision:  2013 July 26

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# int iauDtf2d(const char *scale, int iy, int im, int id,
#              int ihr, int imn, double sec, double *d1, double *d2)
function iauDtf2d(scale::String, iy::Real, im::Real, id::Real, ihr::Real, imn::Real, sec::Real)
   ref_d1 = Ref{Float64}(0.0)
   ref_d2 = Ref{Float64}(0.0)

   status = ccall((:iauDtf2d, libsofa_c), Cint,
                  (Cstring, Cint, Cint, Cint, Cint, Cint, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                  scale, convert(Int32, iy), convert(Int32, im), convert(Int32, id),
                  convert(Int32, ihr), convert(Int32, imn), convert(Float64, sec), ref_d1, ref_d2)

   if status != 0
      @warn "Non-zero return code form iauDtf2d"
   end

   return status, ref_d1[], ref_d2[]
end