text:
module CYCLOPS_v6_2a_PreNPostprocessModule
#export Size_Transform_Seed_Data

export GetEigenGenes
export GetNEigenGenes
export PCA_Transform_Seed_Data
export Repeat_PCA_Transform_Data
export Row_Shuffle
export Bonferroni_Adjust
export Corrected_Cosinor_Statistics
export Corrected_Cosinor_Statistics_Faster,MultiCore_Cosinor_Statistics
export Compile_MultiCore_Cosinor_Statistics
export best_shift_cos, best_shift_cos2
using StatsBase
using MultivariateStats
using Distributions

###################################################
###################################################
function GetEigenGenes(numeric_data::Array{Float64,2},fraction_var=.75, dfrac_var=0.05,maxeig=50)
    u, w, v = svd(numeric_data)
    expvar=cumsum(w.^2)/sum(w.^2)
    ReductionDim1=1+length(expvar[expvar .<= fraction_var]);
    vardif=diff(expvar)
    ReductionDim2=1+length(vardif[vardif .>= dfrac_var]);
    ReductionDim=minimum([ReductionDim1,ReductionDim2])
    ReductionDim=minimum([ReductionDim,maxeig])
    Transform=v[:,1:ReductionDim]'
    ReductionDim, 10*Transform
end

###################################################
###################################################
###################################################
function GetNEigenGenes(numeric_data::Array{Float64,2},nkeep::Integer)
    u, w, v = svd(numeric_data)
    Transform=v[:,1:nkeep]'
    10*Transform
end

###################################################
###################################################
function PCA_Transform_Seed_Data(numeric_data::Array{Float64,2},fraction_var=.75,maxeig=50)
    Mod=fit(PCA,numeric_data;pratio=fraction_var,maxoutdim=maxeig)
    ReductionDim=outdim(Mod)
    Transform=transform(Mod,numeric_data)
    Transform=10*Transform/maximum(Transform)
    ReductionDim, Transform
end

###################################################
###################################################
function Repeat_PCA_Transform_Data(numeric_data::Array{Float64,2},outsize::Integer)
    z_data=zscore(numeric_data,2)
    Mod=fit(PCA,z_data;pratio=1.,maxoutdim=outsize)
    Transform=transform(Mod,z_data)
    Transform
end
##########################################################
###########################################################
function Row_Shuffle(data_matrix::Array{Float64,2})
#ROW SHUFFLE independantly permutes the elements in each row of a matrix
    (nrows,ncols)=size(data_matrix);
    shuffled_matrix=zeros(nrows,ncols)
    for r=1:nrows
        shuffled_matrix[r,:]=data_matrix[r,sample(1:ncols,ncols;replace=false)];
    end
    shuffled_matrix
end

##########################################################
###########################################################
function Bonferroni_Adjust(pvalues::Array{Float64,1})
    n=size(pvalues)[1]
    adjusted=zeros(n)
    for cnt=1:n
        adjusted[cnt]=minimum([n*pvalues[cnt],1.])
    end
    adjusted
end

##########################################################
##########################################################
#####################################################################################################################
function Corrected_Cosinor_Statistics_Faster(expression::Array{Float64,1},o_PREDICTED_PHASELIST::Array{Float64,1},n_shifts=20)

    len=size(o_PREDICTED_PHASELIST)[1]
    shift_list=linspace(2*pi/n_shifts,2*pi,n_shifts)

    min_error=var(expression)*10E20
    best_shift="error"
    for shift in shift_list
        l_PREDICTED_PHASELIST       =mod.((o_PREDICTED_PHASELIST + shift),(2*pi))

        XLINEAR =reshape(l_PREDICTED_PHASELIST,len,1);

        mod_linear          =llsq(XLINEAR,expression);
        m_linear, b_linear  = mod_linear[1:end-1], mod_linear[end];
        predict_linear      =XLINEAR*m_linear + b_linear;

        sse_linear          =sum(abs2.(expression-predict_linear));
        if sse_linear < min_error
            min_error=sse_linear
            best_shift=shift
        end
    end

    l_PREDICTED_PHASELIST=mod.((o_PREDICTED_PHASELIST + best_shift),(2*pi))
    SIN_l_PREDICTED_PHASELIST   =sin.(l_PREDICTED_PHASELIST);
    COS_l_PREDICTED_PHASELIST   =cos.(l_PREDICTED_PHASELIST);

    XLINEAR =reshape(l_PREDICTED_PHASELIST,len,1);
    XFULL   =hcat(l_PREDICTED_PHASELIST,SIN_l_PREDICTED_PHASELIST,COS_l_PREDICTED_PHASELIST);

    mod_linear          =llsq(XLINEAR,expression);
    m_linear, b_linear  = mod_linear[1:end-1], mod_linear[end];

    mod_full            =llsq(XFULL,expression);
    m_full, b_full      = mod_full[1:end-1], mod_full[end];

    predict_linear      =XLINEAR*m_linear + b_linear;
    predict_full        =XFULL*m_full + b_full;

    sse_linear          =sum(abs2.(expression-predict_linear));
    sse_full            =sum(abs2.(expression-predict_full));

    f_metric            =((sse_linear-sse_full)/2)/((sse_full)/(len-5)); #p2=4+1 (+1 =shift) ,p1=2+1  (+1 =shift) ,n=len
    pval= 1-cdf(FDist(2,len-5),f_metric) #Note - null hypothesis this satisifes F distribution, with (p2−p1, n−p2) degrees of freedom.

    ###################
    #Get Phase,Amp,Avg
    ####################
    SIN_o_PREDICTED_PHASELIST   =sin.(o_PREDICTED_PHASELIST);
    COS_o_PREDICTED_PHASELIST   =cos.(o_PREDICTED_PHASELIST);
    XCOSINOR                    =hcat(SIN_o_PREDICTED_PHASELIST,COS_o_PREDICTED_PHASELIST);
    mod_cosinor                 =llsq(XCOSINOR,expression);

    m_cosinor, b_cosinor        = mod_cosinor[1:end-1], mod_cosinor[end];

    predict_cosinor             = XCOSINOR*m_cosinor + b_cosinor;

    sse_cosinor                 =sum(abs2.(expression-predict_cosinor));
    mean_expression             =mean(expression);
    sse_base                    =sum(abs2.(expression-mean_expression));

    r2                          =1-(sse_cosinor/sse_base)

    sinterm, costerm  = mod_cosinor[1], mod_cosinor[2];
    acrophase=atan2(sinterm,costerm);
    amp=sqrt(sinterm^2+costerm^2);
    fitavg=b_cosinor[1]
    avg=mean(expression)
    pval, acrophase, amp, fitavg, avg, r2
end

##########################################################
#####################################################################################################################
#####################################################################################################################
function MultiCore_Cosinor_Statistics(data::Array{Float64,2},PREDICTED_PHASELIST::Array{Float64,1},n_shifts=20)
    ngenes,nsamples=size(data)
    PrbPval=zeros(ngenes);
    PrbPhase=zeros(ngenes);
    PrbAmp=zeros(ngenes);
    PrbFitMean=zeros(ngenes);
    PrbMean=zeros(ngenes);
    PrbR2=zeros(ngenes);

    for row=1:ngenes
        PrbPval[row],PrbPhase[row],PrbAmp[row],PrbFitMean[row],PrbMean[row],PrbR2[row]=Corrected_Cosinor_Statistics_Faster(data[row,:],PREDICTED_PHASELIST,n_shifts)
    end
    PrbPval,PrbPhase,PrbAmp,PrbFitMean,PrbMean,PrbR2
end

##############################################
#####################################################################################################################
function best_shift_cos2(acrophaselist1,acrophaselist2,samplephase1,conversion_flag)

    if conversion_flag=="hours"
        llist2=mod.(acrophaselist2,24)*(pi/12)
    elseif conversion_flag=="radians"
        llist2=mod.(acrophaselist2,2*pi)
    else
        println("FLAG ERROR")
    end

    bestdist=10;
    bestacrophaselist1=acrophaselist1
    bestsamplephase1=samplephase1

    for a in linspace(-pi,pi,192)
        llist1=mod.(acrophaselist1+a,2*pi)
        adist=mean(1 - cos.(float(llist1-llist2)))
        if bestdist>adist
            bestdist=adist
            bestacrophaselist1=llist1
            bestsamplephase1=mod.(samplephase1+a,2*pi)
        end
      end

    for a in linspace(-pi,pi,192)
        llist1=mod.(-(acrophaselist1+a),2*pi)
        adist=mean(1 - cos.(float(llist1-llist2)))
        if bestdist>adist
            bestdist=adist
            bestacrophaselist1=llist1
            bestsamplephase1=mod.(-(samplephase1+a),2*pi)

        end

    end

    bestacrophaselist1,bestsamplephase1
end

##########################################################
function best_shift_cos(list1,list2,conversion_flag)
    nlist=length(list1)

    if conversion_flag=="hours"
        llist2=mod.(list2,24)*(pi/12)
    elseif conversion_flag=="radians"
        llist2=mod.(list2,2*pi)
    else
        println("FLAG ERROR")
    end

    bestdist=10;
    bestlist=zeros(nlist);
    n=0;
    for a in linspace(-pi,pi,192)
        n=n+1
        llist1=mod.(list1+a,2*pi)
        runningscore=0
        for i in 1:nlist
            cosd=(1-cos(llist1[i]-llist2[i]))/2
            runningscore=runningscore+cosd
        end
        adist=runningscore/nlist
        if bestdist>adist
            bestdist=adist
            bestlist=llist1
        end
      end

    for a in linspace(-pi,pi,192)
        n=n+1
        llist1=mod.(-(list1+a),2*pi)
        runningscore=0
        for i in 1:nlist
            cosd=(1-cos(llist1[i]-llist2[i]))/2
            runningscore=runningscore+cosd
        end
        adist=runningscore/nlist

        if bestdist>adist
            bestdist=adist
            bestlist=llist1
        end

    end

    bestlist
end
#####################################################################################################################
function Compile_MultiCore_Cosinor_Statistics(annotated_data::Array{Any,2},PREDICTED_PHASELIST::Array{Float64,1},NumericStartCol=4,n_shifts=20)

alldata_annot=annotated_data[1:end ,1:(NumericStartCol-1)];
alldata_data=Array{Float64}(annotated_data[2:end , NumericStartCol:end]);

ngenes=size(alldata_data)[1]
rowbin=Int(floor(ngenes/5));
tic()
estimated_phaselist=vec(PREDICTED_PHASELIST)
c1=@spawn MultiCore_Cosinor_Statistics(alldata_data[(1:(1*rowbin)),:],estimated_phaselist) ;
c2=@spawn MultiCore_Cosinor_Statistics(alldata_data[((1+1*rowbin):(2*rowbin)),:],estimated_phaselist,n_shifts);
c3=@spawn MultiCore_Cosinor_Statistics(alldata_data[((1+2*rowbin):(3*rowbin)),:],estimated_phaselist,n_shifts);
c4=@spawn MultiCore_Cosinor_Statistics(alldata_data[((1+3*rowbin):(4*rowbin)),:],estimated_phaselist,n_shifts);
c5=@spawn MultiCore_Cosinor_Statistics(alldata_data[((1+4*rowbin):ngenes),:],estimated_phaselist,n_shifts);

PrbPval1,PrbPhase1,PrbAmp1,PrbFitMean1,PrbMean1,PrbRsq1=fetch(c1);
PrbPval2,PrbPhase2,PrbAmp2,PrbFitMean2,PrbMean2,PrbRsq2=fetch(c2);
PrbPval3,PrbPhase3,PrbAmp3,PrbFitMean3,PrbMean3,PrbRsq3=fetch(c3);
PrbPval4,PrbPhase4,PrbAmp4,PrbFitMean4,PrbMean4,PrbRsq4=fetch(c4);
PrbPval5,PrbPhase5,PrbAmp5,PrbFitMean5,PrbMean5,PrbRsq5=fetch(c5);
toc()

PrbPval=vcat(PrbPval1,PrbPval2,PrbPval3,PrbPval4,PrbPval5);
PrbPhase=vcat(PrbPhase1,PrbPhase2,PrbPhase3,PrbPhase4,PrbPhase5);
PrbAmp=vcat(PrbAmp1,PrbAmp2,PrbAmp3,PrbAmp4,PrbAmp5);
PrbMean=vcat(PrbMean1,PrbMean2,PrbMean3,PrbMean4,PrbMean5);
PrbFitMean=vcat(PrbFitMean1,PrbFitMean2,PrbFitMean3,PrbFitMean4,PrbFitMean5);
PrbRsq=vcat(PrbRsq1,PrbRsq2,PrbRsq3,PrbRsq4,PrbRsq5);

PrbPtr= (PrbAmp + PrbFitMean) ./  (PrbFitMean - PrbAmp)
PrbBon=Bonferroni_Adjust(PrbPval);

cosinor_output=hcat(PrbPval,PrbBon,PrbPhase,PrbAmp,PrbFitMean,PrbMean,PrbRsq,PrbPtr)
headers=["pval","bon_pval","phase","amp","fitmean","mean","rsq","ptr"]
headers=reshape(headers,1,length(headers))

cosinor_output=vcat(headers,cosinor_output)
cosinor_output=hcat(alldata_annot,cosinor_output)
cosinor_output
end

end

max_stars_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_issues_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_forks_repo_name: AdithyaVenkateshMohan/CYCLOPS


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""TokenReviewStatus is the result of the token authentication request.

    IoK8sApiAuthenticationV1TokenReviewStatus(;
        audiences=nothing,
        authenticated=nothing,
        error=nothing,
        user=nothing,
    )

    - audiences::Vector{String} : Audiences are audience identifiers chosen by the authenticator that are compatible with both the TokenReview and token. An identifier is any identifier in the intersection of the TokenReviewSpec audiences and the token's audiences. A client of the TokenReview API that sets the spec.audiences field should validate that a compatible audience identifier is returned in the status.audiences field to ensure that the TokenReview server is audience aware. If a TokenReview returns an empty status.audience field where status.authenticated is \"true\", the token is valid against the audience of the Kubernetes API server.
    - authenticated::Bool : Authenticated indicates that the token was associated with a known user.
    - error::String : Error indicates that the token couldn't be checked
    - user::IoK8sApiAuthenticationV1UserInfo : User is the UserInfo associated with the provided token.
"""
mutable struct IoK8sApiAuthenticationV1TokenReviewStatus <: SwaggerModel
    audiences::Any # spec type: Union{ Nothing, Vector{String} } # spec name: audiences
    authenticated::Any # spec type: Union{ Nothing, Bool } # spec name: authenticated
    error::Any # spec type: Union{ Nothing, String } # spec name: error
    user::Any # spec type: Union{ Nothing, IoK8sApiAuthenticationV1UserInfo } # spec name: user

    function IoK8sApiAuthenticationV1TokenReviewStatus(;audiences=nothing, authenticated=nothing, error=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("audiences"), audiences)
        setfield!(o, Symbol("audiences"), audiences)
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("authenticated"), authenticated)
        setfield!(o, Symbol("authenticated"), authenticated)
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        validate_property(IoK8sApiAuthenticationV1TokenReviewStatus, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiAuthenticationV1TokenReviewStatus

const _property_map_IoK8sApiAuthenticationV1TokenReviewStatus = Dict{Symbol,Symbol}(Symbol("audiences")=>Symbol("audiences"), Symbol("authenticated")=>Symbol("authenticated"), Symbol("error")=>Symbol("error"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiAuthenticationV1TokenReviewStatus = Dict{Symbol,String}(Symbol("audiences")=>"Vector{String}", Symbol("authenticated")=>"Bool", Symbol("error")=>"String", Symbol("user")=>"IoK8sApiAuthenticationV1UserInfo")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthenticationV1TokenReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthenticationV1TokenReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1TokenReviewStatus[property_name]

function check_required(o::IoK8sApiAuthenticationV1TokenReviewStatus)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1TokenReviewStatus }, name::Symbol, val)
end

ext: jl


text:
export Subject, SubjectFactory

import Base: show, similar

##

struct SubjectListener{I}
    schedulerinstance :: I
    actor
end

Base.show(io::IO, ::SubjectListener) = print(io, "SubjectListener()")

##

"""
    Subject(::Type{D}; scheduler::H = AsapScheduler())

A Subject is a special type of Observable that allows values to be multicasted to many Observers. Subjects are like EventEmitters.
Every Subject is an Observable and an Actor. You can subscribe to a Subject, and you can call `next!` to feed values as well as `error!` and `complete!`.

Note: By convention, every actor subscribed to a Subject observable is not allowed to throw exceptions during `next!`, `error!` and `complete!` calls.
Doing so would lead to undefined behaviour. Use `safe()` operator to bypass this rule.

See also: [`SubjectFactory`](@ref), [`ReplaySubject`](@ref), [`BehaviorSubject`](@ref), [`safe`](@ref)
"""
mutable struct Subject{D, H, I} <: AbstractSubject{D}
    listeners   :: List{SubjectListener{I}}
    scheduler   :: H
    isactive    :: Bool
    iscompleted :: Bool
    isfailed    :: Bool
    lasterror   :: Any

    Subject{D, H, I}(scheduler::H) where { D, H <: AbstractScheduler, I } = new(List(SubjectListener{I}), scheduler, true, false, false, nothing)
end

function Subject(::Type{D}; scheduler::H = AsapScheduler()) where { D, H <: AbstractScheduler }
    return Subject{D, H, instancetype(D, H)}(scheduler)
end

Base.show(io::IO, ::Subject{D, H}) where { D, H } = print(io, "Subject($D, $H)")

Base.similar(subject::Subject{D, H}) where { D, H } = Subject(D; scheduler = similar(subject.scheduler))

##

isactive(subject::Subject)    = subject.isactive
iscompleted(subject::Subject) = subject.iscompleted
isfailed(subject::Subject)    = subject.isfailed
lasterror(subject::Subject)   = subject.lasterror

setinactive!(subject::Subject)       = subject.isactive    = false
setcompleted!(subject::Subject)      = subject.iscompleted = true
setfailed!(subject::Subject)         = subject.isfailed    = true
setlasterror!(subject::Subject, err) = subject.lasterror   = err

##

function on_next!(subject::Subject{D, H, I}, data::D) where { D, H, I }
    for listener in subject.listeners
        scheduled_next!(listener.actor, data, listener.schedulerinstance)
    end
end

function on_error!(subject::Subject, err)
    if isactive(subject)
        setinactive!(subject)
        setfailed!(subject)
        setlasterror!(subject, err)
        for listener in subject.listeners
            scheduled_error!(listener.actor, err, listener.schedulerinstance)
        end
        empty!(subject.listeners)
    end
end

function on_complete!(subject::Subject)
    if isactive(subject)
        setinactive!(subject)
        setcompleted!(subject)
        for listener in subject.listeners
            scheduled_complete!(listener.actor, listener.schedulerinstance)
        end
        empty!(subject.listeners)
    end
end

##

function on_subscribe!(subject::Subject{D}, actor) where { D }
    if isfailed(subject)
        error!(actor, lasterror(subject))
        return SubjectSubscription(nothing)
    elseif iscompleted(subject)
        complete!(actor)
        return SubjectSubscription(nothing)
    else
        instance = makeinstance(D, subject.scheduler)
        return scheduled_subscription!(subject, actor, instance)
    end
end

function on_subscribe!(subject::Subject, actor, instance)
    listener      = SubjectListener(instance, actor)
    listener_node = pushnode!(subject.listeners, listener)
    return SubjectSubscription(listener_node)
end

##

mutable struct SubjectSubscription <: Teardown
    listener_node :: Union{Nothing, ListNode}
end

as_teardown(::Type{ <: SubjectSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::SubjectSubscription)
    if subscription.listener_node !== nothing
        remove(subscription.listener_node)
        subscription.listener_node = nothing
    end
    return nothing
end

Base.show(io::IO, ::SubjectSubscription) = print(io, "SubjectSubscription()")

##

"""
    SubjectFactory(scheduler::H) where { H <: AbstractScheduler }

A base subject factory that creates an instance of Subject with specified scheduler.

See also: [`AbstractSubjectFactory`](@ref), [`Subject`](@ref)
"""
struct SubjectFactory{ H <: AbstractScheduler } <: AbstractSubjectFactory
    scheduler :: H
end

create_subject(::Type{L}, factory::SubjectFactory) where L = Subject(L, scheduler = similar(factory.scheduler))

Base.show(io::IO, ::SubjectFactory{H}) where H = print(io, "SubjectFactory($H)")


text:
include("get_atmos_ref_states.jl")
using JLD2
using Pkg.Artifacts
using ClimateMachine.ArtifactWrappers
using ClimateMachine.Thermodynamics
const TD = Thermodynamics

@testset "Hydrostatic reference states - regression test" begin
    ref_state_dataset = ArtifactWrapper(
        joinpath(@__DIR__, "Artifacts.toml"),
        "ref_state",
        ArtifactFile[ArtifactFile(
            url = "https://caltech.box.com/shared/static/gyq292ns79wm9xpmy1sse3qtnpcxw54q.jld2",
            filename = "ref_state.jld2",
        ),],
    )
    ref_state_dataset_path = get_data_folder(ref_state_dataset)
    data_file = joinpath(ref_state_dataset_path, "ref_state.jld2")

    RH = 0.5
    (nelem_vert, N_poly) = (20, 4)
    solver_config = get_atmos_ref_states(nelem_vert, N_poly, RH)
    dons_arr = dict_of_nodal_states(solver_config, (Auxiliary(),))
    T = dons_arr["ref_state.T"]
    p = dons_arr["ref_state.p"]
    ρ = dons_arr["ref_state.ρ"]

    @load "$data_file" T_ref p_ref ρ_ref
    @test all(T .≈ T_ref)
    @test all(p .≈ p_ref)
    @test all(ρ .≈ ρ_ref)
end

@testset "Hydrostatic reference states - correctness" begin

    RH = 0.5
    # Fails on (80, 1)
    for (nelem_vert, N_poly) in [(40, 2), (20, 4)]
        solver_config = get_atmos_ref_states(nelem_vert, N_poly, RH)
        dons_arr = dict_of_nodal_states(solver_config)
        phase_type = PhaseEquil
        T = dons_arr["ref_state.T"]
        p = dons_arr["ref_state.p"]
        ρ = dons_arr["ref_state.ρ"]
        q_tot = dons_arr["ref_state.ρq_tot"] ./ ρ
        q_pt = PhasePartition.(q_tot)

        # TODO: test that ρ and p are in discrete hydrostatic balance

        # Test state for thermodynamic consistency (with ideal gas law)
        T_igl =
            TD.air_temperature_from_ideal_gas_law.(Ref(param_set), p, ρ, q_pt)
        @test all(T .≈ T_igl)

        # Test that relative humidity in reference state is approximately
        # input relative humidity
        RH_ref = relative_humidity.(Ref(param_set), T, p, Ref(phase_type), q_pt)
        @show max(abs.(RH .- RH_ref)...)
        @test all(isapprox.(RH, RH_ref, atol = 0.05))
    end

end

max_forks_repo_name: jm-c/CLIMA


text:
"""
	$(SIGNATURES)

Set plot defaults.
"""
function plot_defaults()
    set_theme!(default_theme());
end

"""
	$(SIGNATURES)

Default theme. Note that `Makie` also exports `default_theme`.
"""
function default_theme()
    thm = theme_ggplot2();
    update_theme!(thm;
        # Lines = (linewidth = 3),  # causes errors
        linewidth = 3,
        colormap = default_color_scheme(),
        # fontsize = 10,
        whiskerwidth = 8
        );
    return thm
end

default_color_scheme() = ColorSchemes.corkO;

"""
	$(SIGNATURES)

Color used for a single color graph.
"""
main_color() = get_colors(1, 5);

get_colors(idx) = default_color_scheme()[idx];

"""
	$(SIGNATURES)

Get color `idx` out of `nMax` along the color range defined by the colormap.
"""
function get_colors(idx, nMax :: Integer)
    # nColors = length(default_color_scheme());
    # idxV = round.(Int, LinRange(1, nColors, nMax));
    idxV = 0.05 .+ (0.95 - 0.05) .* idx ./ nMax;
    return get_colors(idxV);
    # return get_colors(idxV[idx])
end

# -----------


text:
@testset "SplayTree" begin
    @testset "pushing values" begin
        t = SplayTree{Int}()
        for i in 1:100
            push!(t, i)
        end

        @test length(t) == 100

        for i in 1:100
            @test haskey(t, i)
        end

        for i = 101:200
            @test !haskey(t, i)
        end
    end

    @testset "deleting values" begin
        t = SplayTree{Int}()
        for i in 1:100
            push!(t, i)
        end
        for i in 1:2:100
            delete!(t, i)
        end

        @test length(t) == 50

        for i in 1:100
            if iseven(i)
                @test haskey(t, i)
            else
                @test !haskey(t, i)
            end
        end

        for i in 1:2:100
            push!(t, i)
        end

        @test length(t) == 100
    end

    @testset "handling different cases of delete!" begin
        t2 = SplayTree()
        for i in 1:100000
            push!(t2, i)
        end

        @test length(t2) == 100000

        nums = rand(1:100000, 8599)
        visited = Set()
        for num in nums
            if !(num in visited)
                delete!(t2, num)
                push!(visited, num)
            end
        end

        for i in visited
            @test !haskey(t2, i)
        end
        @test (length(t2) + length(visited)) == 100000
    end

    @testset "handling different cases of push!" begin
        nums = rand(1:100000, 1000)
        t3 = SplayTree()
        uniq_nums = Set(nums)
        for num in nums
            push!(t3, num)
        end
        @test length(t3) == length(uniq_nums)
    end

    @testset "in" begin
        t4 = SplayTree{Char}()
        push!(t4, 'a')
        push!(t4, 'b')
        @test length(t4) == 2
        @test in('a', t4)
        @test !in('c', t4)
    end

    @testset "search_node" begin
        t5 = SplayTree()
        for i in 1:32
            push!(t5, i)
        end
        n1 = search_node(t5, 21)
        @test n1.data == 21
        n2 = search_node(t5, 35)
        @test n2.data == 32
        n3 = search_node(t5, 0)
        @test n3.data == 1
    end

    @testset "getindex" begin
        t6 = SplayTree{Int}()
        for i in 1:10
            push!(t6, i)
        end
        for i in 1:10
            @test t6[i] == i
        end
        @test_throws KeyError getindex(t6, 0)
        @test_throws KeyError getindex(t6, 11)
    end

    @testset "key conversion in push!" begin
        t7 = SplayTree{Int}()
        push!(t7, Int8(1))
        @test length(t7) == 1
        @test haskey(t7, 1)
    end

    @testset "maximum_node" begin
        t8 = SplayTree()
        @test maximum_node(t8.root) == nothing
        for i in 1:32
            push!(t8, i)
        end
        m1 = maximum_node(t8.root)
        @test m1.data == 32
        node = t8.root
        while node.rightChild != nothing
            m = maximum_node(node.rightChild)
            @test m == m1
            node = node.rightChild
        end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct Point3D{T} <: FieldVector{3,T}
	r :: T
	θ :: T
	ϕ :: T
end
Point3D(x...) = Point3D(promote(x...)...)
struct Point2D{T} <: FieldVector{2,T}
	θ :: T
	ϕ :: T
end
Point3D(x::Tuple) = Point3D(x...)
Point2D(x::Tuple) = Point2D(x...)
Point2D(x::T...) where {T} = throw(ArgumentError("Point2D requires two angular coordinates"))
Point3D(x::T...) where {T} = throw(ArgumentError("Point2D requires three angular coordinates"))
Point2D(x...) = Point2D(promote(x...)...)

Point3D(r, p::Point2D) = Point3D(r, p.θ, p.ϕ)
Point3D(p::Point2D) = Point3D(1, p)
Point2D(p::Point3D) = Point2D(p.θ, p.ϕ)
Point2D(p::Point2D) = p
Point3D(p::Point3D) = p

radiusvector(θ, ϕ) = SVector{3}(sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ))
radiusvector(r, θ, ϕ) = SVector{3}(sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ)) .* r
radiusvector(p::Point3D) = radiusvector(p.r, p.θ, p.ϕ)
radiusvector(p::Point2D) = radiusvector(p.θ, p.ϕ)

Point3D(x::AbstractVector) = Point3D(SVector{3}(x))
Point2D(x::AbstractVector) = Point2D(SVector{3}(x))
function polarangles(v)
    x, y, z = v
    θ = acos(z)
    if abs(z) == 1
        ϕ = 0.0
    else
        ϕ = atan(x, y)
    end
    promote(θ, ϕ)
end
Point3D(v::SVector{3}) = Point3D(norm(v), polarangles(v)...)
Point2D(v::SVector{3}) = Point2D(polarangles(normalize(v))...)

const SphericalPoint = Union{Point2D,Point3D}

Base.broadcastable(p::SphericalPoint) = Ref(p)

for f in [:cosχ, :∂θ₁cosχ, :∂θ₂cosχ, :∂ϕ₁cosχ, :∂ϕ₂cosχ]
    @eval $f(n1, n2) = $f(Point2D(n1), Point2D(n2))
end
cosχ(n1::SphericalPoint, n2::SphericalPoint) = radiusvector(Point2D(n1)) ⋅ radiusvector(Point2D(n2))
∂θ₁cosχ(n1::SphericalPoint, n2::SphericalPoint) = epsilon(cosχ(Point2D(Dual(n1.θ,1), n1.ϕ), n2))
∂θ₂cosχ(n1::SphericalPoint, n2::SphericalPoint) = epsilon(cosχ(n1, Point2D(Dual(n2.θ,1), n2.ϕ)))
∂ϕ₁cosχ(n1::SphericalPoint, n2::SphericalPoint) = epsilon(cosχ(Point2D(n1.θ, Dual(n1.ϕ,1)), n2))
∂ϕ₂cosχ(n1::SphericalPoint, n2::SphericalPoint) = epsilon(cosχ(n1, Point2D(n2.θ, Dual(n2.ϕ, 1))))

shiftϕ(p::Point2D, Δϕ) = Point2D(p.θ, p.ϕ + Δϕ)
shiftϕ(p::Point3D, Δϕ) = Point3D(p.r, p.θ, p.ϕ + Δϕ)


text:
using Test
using Suppressor
using Pkg
using LibGit2
using PackageStates

const SIGNATURE = LibGit2.Signature("DUMMY", "DUMMY@DUMMY.DOM", round(time()), 0)

function mkdummypackage(tmpdir, name)
    targetdir = joinpath(tmpdir, name)
    cp(joinpath(@__DIR__, "test_packages", name), targetdir)
    LibGit2.with(LibGit2.init(targetdir)) do repo
        LibGit2.add!(repo, "*")
        LibGit2.commit(repo, "msg"; author=SIGNATURE, committer=SIGNATURE)
    end
    return targetdir
end

remove_dateline_from_diff(diffstr) = join(split(diffstr, "\n")[union(1:4,6:end)], "\n")
remove_dateline_and_header_from_diff(diffstr) = join(split(diffstr, "\n")[union(1:2,4,6:end)], "\n")

@testset "Basics" begin
    mktempdir() do tmp

        env1 = mkdir(joinpath(tmp, "env1"))
        env2 = mkdir(joinpath(tmp, "env2"))

        dummy = mkdummypackage(tmp, "DummyPackage")
        dummy2 = mkdummypackage(tmp, "DummyPackage_v2")

        th1 = Sys.iswindows() ? "1bd2b16b793dfbf96aef17385f635729ae32a43c" : "dd574217160ae714ff496b9239a5ae1a4d819aa8"
        th2 = Sys.iswindows() ? "90ca0b300186580cfe89f5336ec58453257a1ec6" : "98cea5b18356123cda026692eafb1e4a55813dac"

        Pkg.activate(env1)
        Pkg.add(path=dummy)
        Pkg.activate(env2)
        Pkg.add(path=dummy2)

        Pkg.activate(env1)
        @eval using DummyPackage

        s = state(DummyPackage)
        @test s.id == Base.PkgId(Base.UUID("a5a70863-7a97-4c01-857e-744174fcb92d"), "DummyPackage")
        @test s.project == joinpath(env1, "Project.toml")
        @test s.load_path[1] == joinpath(env1, "Project.toml")
        @test s.tree_hash == th1
        @test s.manifest_tree_hash == s.tree_hash

        Pkg.activate(env2)

        s2 = state(DummyPackage)
        @test s2.tree_hash == th1
        @test s2.manifest_tree_hash == th2
        @test s2.load_path[1] == joinpath(env2, "Project.toml")

        @test @capture_out(diff_states_all(:on_load => :newest)) == ""
        snew = state(DummyPackage, :newest)
        @test snew.manifest_tree_hash == th1

        @test diff_states(DummyPackage, print = false)
        @test @capture_out(diff_states(DummyPackage)) ≠ ""
        @test remove_dateline_from_diff(@capture_out(diff_states(DummyPackage))) == remove_dateline_from_diff(@capture_out(diff_states(DummyPackage, :newest => :current)))

        diff_states_all(print = false, update = true)
        supdated = state(DummyPackage, :newest)
        @test supdated.manifest_tree_hash == th2
        @test @capture_out(diff_states_all(:on_load => :newest)) ≠ ""
        @test remove_dateline_and_header_from_diff(@capture_out(diff_states(DummyPackage, :on_load => :newest))) == remove_dateline_and_header_from_diff(@capture_out(diff_states(DummyPackage, :on_load => :current)))
        Pkg.activate(env1)
        @test diff_states_all(:on_load => :current, print=false) == [PackageStates]

        @test recorded_modules() == Set([DummyPackage, PackageStates])
    end
end


text:
### abstract type hierarchy

abstract type Basis end
abstract type Basislike <: Basis end

struct DerivBasis <: Basislike end # Basis in which derivatives are sparse for a given field
const Ð! = DerivBasis
const Ð  = DerivBasis
struct LenseBasis <: Basislike end # Basis in which lensing is a pixel remapping for a given field
const Ł! = LenseBasis
const Ł  = LenseBasis

## fields
abstract type Field{B<:Basis,T} <: AbstractVector{T} end
abstract type ImplicitField{B<:Basis,T} <: Field{B,T} end

## linear operators
abstract type ImplicitOp{T} <: AbstractMatrix{T} end
const DiagOp{F<:Field,T} = Diagonal{T,F}
# includes Diagonal and Adjoint wrappers, which can't be made <:ImplicitOp directly
const FieldOp{T} = Union{ImplicitOp{T},Adjoint{T,<:ImplicitOp{T}},Diagonal{T,<:Field{<:Any,T}}}

# scalars which are allowed in our expressions must be real because we
# implicitly assume our maps are real, and addition/multiplication by a complex
# number, even of the fourier transform, would break this.
const Scalar = Real

# other useful unions
const FieldOrOp   = Union{Field,FieldOp}
const FieldOpScal = Union{Field,FieldOp,Scalar}

## basis types
abstract type BasisProd{Bs} <: Basis end

struct Basis2Prod{B₁,B₂}    <: BasisProd{Tuple{B₁,B₂}} end
struct Basis3Prod{B₁,B₂,B₃} <: BasisProd{Tuple{B₁,B₂,B₃}} end

struct Map     <: Basis end
struct Fourier <: Basis end
struct 𝐈       <: Basis end
struct 𝐐𝐔      <: Basis end
struct 𝐄𝐁      <: Basis end

const QUMap      = Basis2Prod{    𝐐𝐔, Map     }
const EBMap      = Basis2Prod{    𝐄𝐁, Map     }
const QUFourier  = Basis2Prod{    𝐐𝐔, Fourier }
const EBFourier  = Basis2Prod{    𝐄𝐁, Fourier }
const IQUMap     = Basis3Prod{ 𝐈, 𝐐𝐔, Map     }
const IEBMap     = Basis3Prod{ 𝐈, 𝐄𝐁, Map     }
const IQUFourier = Basis3Prod{ 𝐈, 𝐐𝐔, Fourier }
const IEBFourier = Basis3Prod{ 𝐈, 𝐄𝐁, Fourier }

# handy for picking out anything Map/Fourier
const SpatialBasis{B,I,P} = Union{B, Basis2Prod{P,B}, Basis3Prod{I,P,B}}

# for printing
for F in ["QUMap", "EBMap", "QUFourier", "EBFourier", "IQUMap", "IEBMap", "IQUFourier", "IEBFourier"]
    @eval typealias(::Type{$(Symbol(F))}) = $F
end

### basis

## generic promotion rules which might change basis
basis_promotion_rules = Dict(
    # spin-0
    (Map,        Fourier)     => Map,
    # spin-2
    (QUMap,      QUFourier)   => QUMap,
    (EBMap,      EBFourier)   => EBFourier,
    (QUMap,      EBMap)       => QUMap,
    (QUFourier,  EBFourier)   => QUFourier,
    (QUMap,      EBFourier)   => QUMap,
    (QUFourier,  EBMap)       => QUFourier,
    # spin-(0,2)
    (IQUMap,     IQUFourier)  => IQUMap,
    (IEBMap,     IEBFourier)  => IEBFourier,
    (IQUMap,     IEBMap)      => IQUMap,
    (IQUFourier, IEBFourier)  => IQUFourier,
    (IQUMap,     IEBFourier)  => IQUMap,
    (IQUFourier, IEBMap)      => IQUFourier
)
for ((B₁,B₂),B) in basis_promotion_rules
    @eval promote_basis_generic_rule(::$B₁, ::$B₂) = $B()
end
promote_basis_generic_rule(::B, ::B) where {B<:Basis} = B()
promote_basis_generic_rule(::Any, ::Any) = Unknown()
promote_basis_generic(x, y) = select_known_rule(promote_basis_generic_rule, x, y)
unknown_rule_error(::typeof(promote_basis_generic_rule), ::B₁, ::B₂) where {B₁, B₂} =
    error("Can't promote fields in $(typealias(B₁)) and $(typealias(B₂)) bases.")

## stricter rules used only in broadcasting
promote_basis_strict_rule(::B,   ::B )        where {B <:Basis}                                                    = B()
promote_basis_strict_rule(::B,   ::Basislike) where {B <:Basis}                                                    = B()
promote_basis_strict_rule(::B₂,  ::B₀)        where {B₀<:Union{Map,Fourier}, B₂ <:Basis2Prod{  <:Union{𝐐𝐔,𝐄𝐁},B₀}} = B₂()
promote_basis_strict_rule(::B₀₂, ::B₀)        where {B₀<:Union{Map,Fourier}, B₀₂<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁},B₀}} = B₀₂()
promote_basis_strict_rule(::Any, ::Any)                                                                            = Unknown()
promote_basis_strict(x, y) = select_known_rule(promote_basis_strict_rule, x, y)
unknown_rule_error(::typeof(promote_basis_strict_rule), ::B₁, ::B₂) where {B₁, B₂} =
    error("Can't broadcast fields in $(typealias(B₁)) and $(typealias(B₂)) bases.")

## applying bases

# Map(B) or Fourier(B) for another basis, B
(::Type{B})(::Type{B′}) where {B<:Union{Map,Fourier}, B′<:Union{Map,Fourier}} = B
(::Type{B})(::Type{Basis2Prod{  Pol,B′}}) where {Pol, B<:Union{Map,Fourier}, B′<:Union{Map,Fourier}} = Basis2Prod{  Pol,B}
(::Type{B})(::Type{Basis3Prod{𝐈,Pol,B′}}) where {Pol, B<:Union{Map,Fourier}, B′<:Union{Map,Fourier}} = Basis3Prod{𝐈,Pol,B}

# A "basis-like" object, e.g. the lensing basis Ł or derivative basis Ð. For any
# particular types of fields, these might be different actual bases, e.g. the
# lensing basis is Map for S0 but QUMap for S2.
(::Type{B})(f::F) where {F<:Field,B<:Basislike} = B(F)(f)
(::Type{B})(f′::Field, f::F) where {F<:Field,B<:Basislike} = B(F)(f′,f)
(::Type{B})(a::AbstractArray{<:Field}...) where {B<:Basis} = B.(a...)

# B(f) where B is a basis converts f to that basis. This is the fallback if the
# field is already in the right basis.
(::Type{B})(f::Field{B}) where {B<:Basis} = f

# The abstract `Basis` type means "any basis", hence this conversion rule:
Basis(f::Field) = f

# B(f′, f) converts f to basis B and stores the result inplace in f′. If f is
# already in basis B, we just return f (but note, we never actually set f′ in
# this case, which is more efficient, but necessitates some care when using this
# construct)
(::Type{B})(f′::Field{B}, f::Field{B}) where {B<:Basis} = f

# convenience "getter" functions for the B basis type parameter
basis(f::F) where {F<:Field} = basis(F)
basis(::Type{<:Field{B}}) where {B<:Basis} = B
basis(::Type{<:Field}) = Basis
basis(::AbstractVector) = Basis

### printing
typealias(::Type{B}) where {B<:Basis} = string(B.name.name)
Base.show_datatype(io::IO, t::Type{<:Union{Field,FieldOp}}) = print(io, typealias(t))
Base.isempty(::ImplicitOp) = true
Base.isempty(::ImplicitField) = true
Base.isempty(::Diagonal{<:Any,<:ImplicitField}) = true
function Base.summary(io::IO, x::FieldOp)
    try
        print(io, join(size(x), "×"), " ")
    catch err
        @assert err isa MethodError
        print(io, "⍰×⍰ ")
    end
    Base.showarg(io, x, true)
end
function Base.summary(io::IO, x::ImplicitField)
    try
        print(io, size(x,1), "-element ")
    catch err
        @assert err isa MethodError
        print(io, "⍰-element ")
    end
    Base.showarg(io, x, true)
end

# without this, *sometimes* IJulia doesnt print the field types right, but I dont really understand it
@init @require IJulia="7073ff75-c697-5162-941a-fcdaad2a7d2a" begin
    Base.show(io::IOContext{IOBuffer}, t::Type{<:Union{Field,FieldOp}}) = print(io, typealias(t))
end

### logdet

"""
    logdet(L::FieldOp, θ)

If `L` depends on `θ`, evaluates `logdet(L(θ))` offset by its fiducial value at
`L()`. Otherwise, returns 0.
"""
logdet(L::FieldOp, θ) = depends_on(L,θ) ? logdet(L()\L(θ)) : 0
logdet(L::Int, θ) = 0 # the default returns Float64 which unwantedly poisons the backprop to Float64
logdet(L, θ) = logdet(L)

### Simulation

@doc doc"""
    simulate(Σ; rng=global_rng_for(Σ), seed=nothing)

Draw a simulation from the covariance matrix `Σ`, i.e. draw a random vector
$\xi$ such that the covariance $\langle \xi \xi^\dagger \rangle = \Sigma$.

The random number generator `rng` will be used and advanced in the proccess, and
is by default the appropriate one depending on if `Σ` is backed by `Array` or
`CuArray`.

The `seed` argument can also be used to seed the `rng`.
"""
function simulate(Σ; rng=global_rng_for(Σ), seed=nothing, kwargs...)
    isnothing(seed) || seed!(rng, seed)
    simulate(rng, Σ; kwargs...)
end
function white_noise(ξ; rng=global_rng_for(ξ), seed=nothing, kwargs...)
    isnothing(seed) || seed!(rng, seed)
    white_noise(ξ, rng; kwargs...)
end
function fixed_white_noise(ξ; rng=global_rng_for(Σ), seed=nothing, kwargs...)
    isnothing(seed) || seed!(rng, seed)
    fixed_white_noise(ξ, rng; kwargs...)
end

global_rng_for(x::T) where {T<:AbstractArray} = global_rng_for(T)
global_rng_for(::Type{<:Array}) = Random.default_rng()

### spin adjoints
# represents a field which is adjoint over just the "spin" indices.
# multiplying such a field by a non-adjointed one should be the inner
# product over just the spin indices, hence return a spin-0 field,
# rather than a scalar. note: these are really only lightly used in
# one place in LenseFlow, so they have almost no real functionality
# for now.
struct SpinAdjoint{F<:Field}
    f :: F
end
spin_adjoint(f::Field) = SpinAdjoint(f)

### misc

function ≈(x::Field{<:Any,T}, y::Field{<:Any,S}; atol=0, rtol=Base.rtoldefault(T,S,atol)) where {T,S}
    _norm(x) = norm(unbatch(norm(x)))
    _norm(x - y) <= max(atol, rtol*max(_norm(x), _norm(y)))
end

# allows L(θ) or L(ϕ) to work and be a no-op for things which don't
# depend on parameters or a field
(L::Union{FieldOp,UniformScaling})(::Union{Field,NamedTuple}) = L

# allow using `I` as a lensing operator to represent no lensing
alloc_cache(x, ::Any) = x
cache(x, ::Any) = x
cache!(x, ::Any) = x

# caching for adjoints
cache(L::Adjoint, f) = cache(L',f)'
cache!(L::Adjoint, f) = cache!(L',f)'

# todo: fix this
*(::UniformScaling{Bool}, L::FieldOp) = L

# we use Field cat'ing mainly for plotting, e.g. plot([f f; f f]) plots a 2×2
# matrix of maps. the following definitions make it so that Fields aren't
# splatted into a giant matrix when doing [f f; f f] (which they would othewise
# be since they're Arrays)
hvcat(rows::Tuple{Vararg{Int}}, values::Field...) = hvcat(rows, ([x] for x in values)...)
hcat(values::Field...) = hcat(([x] for x in values)...)

### printing
print_array(io::IO, f::Field) = !isempty(f) && print_array(io, f[:])
show_vector(io::IO, f::Field) = !isempty(f) && show_vector(io, f[:])
Base.has_offset_axes(::Field) = false # needed for Diagonal(::Field) if the Field is implicitly-sized

# addition/subtraction works between any fields and scalars, promotion is done
# automatically if fields are in different bases
for op in (:+,:-), (T1,T2,promote) in ((:Field,:Scalar,false),(:Scalar,:Field,false),(:Field,:Field,true))
    @eval ($op)(a::$T1, b::$T2) = broadcast($op, ($promote ? promote(a,b) : (a,b))...)
end

# multiplication/division is not strictly defined for abstract vectors, but
# make it work anyway if the two fields are exactly the same type, in which case
# its clear we wanted broadcasted multiplication/division.
for op in (:*, :/)
    @eval ($op)(a::Field{B}, b::Field{B}) where {B} = broadcast($op, a, b)
end

# needed unless I convince them to undo the changes here:
# https://github.com/JuliaLang/julia/pull/35257#issuecomment-657901503
if VERSION>v"1.4"
    *(x::Adjoint{<:Number,<:Field}, y::Field) = dot(x.parent,y)
end

# package code should use batch_index(f, 1), but for interactive work its
# convenient to be able to do f[!,1]
getindex(x::Union{Real,Field,FieldOp}, ::typeof(!), I) = batch_index(x, I)

# adapt_structure(to, x::Union{ImplicitOp,ImplicitField}) = x
# adapt_structure(to, arr::Array{<:Field}) = map(f -> adapt(to,f), arr)

one(f::Field) = fill!(similar(f), one(eltype(f)))
norm(f::Field) = sqrt(dot(f,f))
# sum_kbn(f::Field) = sum_kbn(f[:])

@init @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
    # never try to auto-convert Fields or FieldOps to Python arrays
    PyCall.PyObject(x::Union{FieldOp,Field}) = PyCall.pyjlwrap_new(x)
end


text:
# https://fluxml.ai/Zygote.jl/dev/adjoints/#Gradient-Reflection-1

# Treat everything immutable during differentiation:
@adjoint possible(_args...) = false, _ -> nothing


text:
export add_find_max, add_find_min, add_find_ith

function add_find_max(mgr::Ptr{Manager}, f::Ptr{Node})
    result = ccall((:Cudd_addFindMax, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}), mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_find_min(mgr::Ptr{Manager}, f::Ptr{Node})
    result = ccall((:Cudd_addFindMin, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}), mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_find_ith(mgr::Ptr{Manager}, f::Ptr{Node}, i::Int)
    result = ccall((:Cudd_addIthBit, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Cuint), mgr, f, i)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
module Sigil

# using DifferentialEquations
# using Plots

mutable struct Star
    M::AbstractFloat
    m::AbstractFloat
    mf::AbstractFloat
    X::AbstractFloat
    Y::AbstractFloat
    grid::Array{AbstractFloat}
    r::Array{AbstractFloat}
    l::Array{AbstractFloat}
    P::Array{AbstractFloat}
    T::Array{AbstractFloat}
end

Star(M) = Star(M, 0.01 * Msun, 0.8 * M, 0.7, 0.28,
               zeros(1), zeros(1), zeros(1), zeros(1), zeros(1))
Star(M, X, Y) = Star(M, 0.01 * Msun, 0.8 * M, X, Y,
                     zeros(1), zeros(1), zeros(1), zeros(1), zeros(1))

export score,
    init_guess,
    deriv!,
    Msun,
    Rsun,
    Star,
    load1,
    load2,
    shootf!,
    plot_star

include("constants.jl")
include("opacity.jl")
include("equations.jl")
include("physics.jl")
include("shootf.jl")
include("plot.jl")

end


text:
using SIMDPoly
using BenchmarkTools
using PyPlot

NSAMPLES = 100_000

degrees = 3:40
simd_times, ref_times = [], []

for n in degrees
    poly1 = ntuple(x -> Float64(x), n)
    packed = packpoly1x4r(poly1)

    our = @benchmark evalpoly1x4r(z, $packed) samples=NSAMPLES setup=(z=rand())
    ref = @benchmark evalpoly(z, $poly1) samples=NSAMPLES setup=(z=rand())
    push!(simd_times, minimum(our).time)
    push!(ref_times, minimum(ref).time)
    #println(minimum(our).time, " / ", minimum(ref).time)
end

ioff()
plot(degrees, ref_times, linewidth=3, label="reference (Base.Math.evalpoly)")
plot(degrees, simd_times, linewidth=3, label="SIMD (4 x f64 lanes)")
legend()
ylim(ymin=0)
grid(axis="y", alpha=0.5)
xlabel("N")
ylabel("time (ns)")
title("degree N polynomial evaluation, real input")
savefig("bm_plot_real.svg")


text:
function clone(obj::NumberFormat)
    return jcall(obj, "clone", Object, ())
end

function equals(obj::NumberFormat, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function format(obj::AbstractFormat, arg0::jdouble, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (jdouble, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::AbstractFormat, arg0::jlong, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (jlong, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::Format, arg0::Object)
    return jcall(obj, "format", JString, (Object,), arg0)
end

function format(obj::NumberFormat, arg0::Object, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (Object, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::NumberFormat, arg0::jdouble)
    return jcall(obj, "format", JString, (jdouble,), arg0)
end

function format(obj::NumberFormat, arg0::jlong)
    return jcall(obj, "format", JString, (jlong,), arg0)
end

function format_to_character_iterator(obj::Format, arg0::Object)
    return jcall(obj, "formatToCharacterIterator", AttributedCharacterIterator, (Object,), arg0)
end

function get_available_locales(::Type{NumberFormat})
    return jcall(NumberFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_compact_number_instance(::Type{NumberFormat})
    return jcall(NumberFormat, "getCompactNumberInstance", NumberFormat, ())
end

function get_compact_number_instance(::Type{NumberFormat}, arg0::Locale, arg1::NumberFormat_Style)
    return jcall(NumberFormat, "getCompactNumberInstance", NumberFormat, (Locale, NumberFormat_Style), arg0, arg1)
end

function get_currency(obj::NumberFormat)
    return jcall(obj, "getCurrency", Currency, ())
end

function get_currency_instance(::Type{NumberFormat})
    return jcall(NumberFormat, "getCurrencyInstance", NumberFormat, ())
end

function get_currency_instance(::Type{NumberFormat}, arg0::Locale)
    return jcall(NumberFormat, "getCurrencyInstance", NumberFormat, (Locale,), arg0)
end

function get_denominator_format(obj::AbstractFormat)
    return jcall(obj, "getDenominatorFormat", NumberFormat, ())
end

function get_instance(::Type{NumberFormat})
    return jcall(NumberFormat, "getInstance", NumberFormat, ())
end

function get_instance(::Type{NumberFormat}, arg0::Locale)
    return jcall(NumberFormat, "getInstance", NumberFormat, (Locale,), arg0)
end

function get_integer_instance(::Type{NumberFormat})
    return jcall(NumberFormat, "getIntegerInstance", NumberFormat, ())
end

function get_integer_instance(::Type{NumberFormat}, arg0::Locale)
    return jcall(NumberFormat, "getIntegerInstance", NumberFormat, (Locale,), arg0)
end

function get_maximum_fraction_digits(obj::NumberFormat)
    return jcall(obj, "getMaximumFractionDigits", jint, ())
end

function get_maximum_integer_digits(obj::NumberFormat)
    return jcall(obj, "getMaximumIntegerDigits", jint, ())
end

function get_minimum_fraction_digits(obj::NumberFormat)
    return jcall(obj, "getMinimumFractionDigits", jint, ())
end

function get_minimum_integer_digits(obj::NumberFormat)
    return jcall(obj, "getMinimumIntegerDigits", jint, ())
end

function get_number_instance(::Type{NumberFormat})
    return jcall(NumberFormat, "getNumberInstance", NumberFormat, ())
end

function get_number_instance(::Type{NumberFormat}, arg0::Locale)
    return jcall(NumberFormat, "getNumberInstance", NumberFormat, (Locale,), arg0)
end

function get_numerator_format(obj::AbstractFormat)
    return jcall(obj, "getNumeratorFormat", NumberFormat, ())
end

function get_percent_instance(::Type{NumberFormat})
    return jcall(NumberFormat, "getPercentInstance", NumberFormat, ())
end

function get_percent_instance(::Type{NumberFormat}, arg0::Locale)
    return jcall(NumberFormat, "getPercentInstance", NumberFormat, (Locale,), arg0)
end

function get_rounding_mode(obj::NumberFormat)
    return jcall(obj, "getRoundingMode", RoundingMode, ())
end

function hash_code(obj::NumberFormat)
    return jcall(obj, "hashCode", jint, ())
end

function is_grouping_used(obj::NumberFormat)
    return jcall(obj, "isGroupingUsed", jboolean, ())
end

function is_parse_integer_only(obj::NumberFormat)
    return jcall(obj, "isParseIntegerOnly", jboolean, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function parse(obj::NumberFormat, arg0::JString)
    return jcall(obj, "parse", Number, (JString,), arg0)
end

function parse(obj::NumberFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", Number, (JString, ParsePosition), arg0, arg1)
end

function parse_object(obj::Format, arg0::JString)
    return jcall(obj, "parseObject", Object, (JString,), arg0)
end

function parse_object(obj::NumberFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parseObject", Object, (JString, ParsePosition), arg0, arg1)
end

function set_currency(obj::NumberFormat, arg0::Currency)
    return jcall(obj, "setCurrency", void, (Currency,), arg0)
end

function set_grouping_used(obj::NumberFormat, arg0::jboolean)
    return jcall(obj, "setGroupingUsed", void, (jboolean,), arg0)
end

function set_maximum_fraction_digits(obj::NumberFormat, arg0::jint)
    return jcall(obj, "setMaximumFractionDigits", void, (jint,), arg0)
end

function set_maximum_integer_digits(obj::NumberFormat, arg0::jint)
    return jcall(obj, "setMaximumIntegerDigits", void, (jint,), arg0)
end

function set_minimum_fraction_digits(obj::NumberFormat, arg0::jint)
    return jcall(obj, "setMinimumFractionDigits", void, (jint,), arg0)
end

function set_minimum_integer_digits(obj::NumberFormat, arg0::jint)
    return jcall(obj, "setMinimumIntegerDigits", void, (jint,), arg0)
end

function set_parse_integer_only(obj::NumberFormat, arg0::jboolean)
    return jcall(obj, "setParseIntegerOnly", void, (jboolean,), arg0)
end

function set_rounding_mode(obj::NumberFormat, arg0::RoundingMode)
    return jcall(obj, "setRoundingMode", void, (RoundingMode,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using IterativeClosestPoint
using CoordinateTransformations
using Rotations
using Distances
using StaticArrays
using LinearAlgebra
using Random
using Test

function randomtform(; s=10, maxangle=π)
    randomaxis = normalize(rand(3) .- 0.5)
    randomangle = (rand() .- 0.5) * 2 * maxangle
    randomrotation = AngleAxis(randomangle, randomaxis...)
    randomtranslation = s * (rand(3) .- 0.5)
    return AffineMap(randomrotation, randomtranslation)
end

@testset "Kabsch" begin
    for i=1:10
        x = SMatrix{3,20}(20 * (rand(3, 20) .- 0.5))
        randtform = randomtform()
        y = randtform(x)

        tform = kabsch(x, y)
        @test tform.linear ≈ randtform.linear
        @test tform.translation ≈ randtform.translation
    end
end

@testset "point-to-point ICP" begin
     # pure shuffling of points (no transformation)
    for i=1:10
        x = SMatrix{3,20}(20 * (rand(3, 20) .- 0.5))
        perm = randperm(20)
        y = (hcat([x[:,i] for i in perm]...))

        matches = iterate_kabsch(y,x)
        @test perm == [m[2] for m in matches]
    end

    # small-rotation with translation
    for i=1:10
        x = SMatrix{3,20}(20 * (rand(3, 20) .- 0.5))
        randtform = randomtform(maxangle = π/4)
        y = randtform(x)
        perm = randperm(20)
        y = (hcat([y[:,i] for i in perm]...))

        matches = iterate_kabsch(y,x)
        @test perm == [m[2] for m in matches]
    end
end


text:
export move, normal, normfrom, direction, vector, extend

xs(l::Line) = SVector(l.from.x, l.to.x)
ys(l::Line) = SVector(l.from.y, l.to.y)

function intersection(l1::Line, l2::Line)
    x1, x2 = xs(l1)
    y1, y2 = ys(l1)
    x3, x4 = xs(l2)
    y3, y4 = ys(l2)
    t = ((x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)) /
        ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))

    px = x1 + t * (x2 - x1)
    py = y1 + t * (y2 - y1)
    Point(px, py)
end

vector(l::Line) = from_to(l.from, l.to)
Base.angle(l::Line) = angle(vector(l))
Base.length(l::Line) = magnitude(vector(l))
fraction(l::Line, frac::Real) = between(l.from, l.to, frac)
center(l::Line) = fraction(l, 0.5)
reversed(l::Line) = Line(l.to, l.from)
direction(l::Line) = normalize(vector(l))
normal(l::Line, reverse::Bool=false) = rotate(direction(l), reverse ? deg(-90) : deg(90))
normal(l::Line, length::Real) = normal(l, false) * length
normfrom(l::Line, frac::Real, length::Real) = fraction(l, frac) + normal(l, length)
move(l::Line, p::Point) = Line(l.from + p, l.to + p)
Base.:+(l::Line, p::Point) = move(l, p)
Base.:+(p::Point, l::Line) = move(l, p)
Base.:-(l::Line, p::Point) = move(l, -p)

Base.:*(l::Line, r::Real) = Line(r * l.from, r * l.to)
Base.:/(l::Line, r::Real) = Line(l.from / r, l.to / r)

function scale(l::Line, scalar::Real)
    movement = (scalar - 1) * vector(l) / 2
    Line(l.from - movement, l.to + movement)
end

function scaleto(l::Line, len::Real)
    scalar = len / length(l)
    scale(l, scalar)
end

function extend(l::Line, len::Real, fraction_forward::Real=1)
    vec = len * direction(l)
    movement_to = fraction_forward * vec
    movement_from = (1 - fraction_forward) * vec
    Line(l.from - movement_from, l.to + movement_to)
end

function rotate(l::Line, angle::Angle; around::Point=Point(0, 0))
    Line(
        rotate(l.from, angle, around=around),
        rotate(l.to, angle, around=around)
    )
end

needed_attributes(::Type{Line}) = (Visible, Linewidth, Stroke, Linestyle)

# needed_attributes(::Type{LineSegments}) = (Visible, Linewidth, Stroke, Linestyle)
#
# Base.convert(::Type{LineSegments}, ls::Vector{Line}) = LineSegments(ls)

function bbox(l::Line)
    mi = min.(l.from.xy, l.to.xy)
    ma = max.(l.from.xy, l.to.xy)
    BBox(P(mi...), P(ma...))
end


text:
using MultivariateFunctions
using DataStructures: OrderedDict
using DataFrames

# Chebyshev approximation
func = sin
nodes  =  12
degree =  8
left   = -2.0
right  =  5.0
limits =  OrderedDict{Symbol,Tuple{Float64,Float64}}([:default] .=> [(left, right)])
approxim = create_chebyshev_approximation(func, nodes, degree, limits)
X = convert(Array{Float64,1}, left:0.01:right)
y = func.(X)
y_approx = evaluate.(Ref(approxim), X)
maximum(abs.(y .- y_approx)) < 0.01

func = exp
nodes  =  12
degree =  8
left   =  1.0
right  =  5.0
limits =  OrderedDict{Symbol,Tuple{Float64,Float64}}([:default] .=> [(left, right)])
approxim = create_chebyshev_approximation(func, nodes, degree, limits)
X = convert(Array{Float64,1}, left:0.01:right)
y = func.(X)
y_approx = evaluate.(Ref(approxim), X)
maximum(abs.(y .- y_approx)) < 0.01

function func1(a::Float64,b::Float64,c::Float64)
    return sin(a)* c + a * sqrt(b)
end
func = func1
nodes  =  8
degree =  4
function_takes_Dict = false
limits = OrderedDict{Symbol,Tuple{Float64,Float64}}([:a, :b, :c] .=> [(-2.0,0.5), (0.0,4.0), (5.0,11.0)])
approxim = create_chebyshev_approximation(func, nodes, degree, limits, function_takes_Dict)
st =  0.99
d = vcat(Iterators.product(convert(Array{Float64,1}, limits[:a][1]:st:limits[:a][2]),convert(Array{Float64,1}, limits[:b][1]:st:limits[:b][2]),convert(Array{Float64,1}, limits[:c][1]:st:limits[:c][2]))...)
dd = DataFrame()
dd[:a], dd[:b], dd[:c]  = vcat.(d...)
y = Array{Float64,1}(undef, size(dd)[1])
for i in 1:length(y)
    y[i] = func.(dd[i,:a], dd[i,:b], dd[i,:c])
end
y_approx = evaluate(approxim, dd)
maximum(abs.(y .- y_approx)) < 0.5
# This is commented out to avoid the extra dependency. It should be true though.
#using Distributions
#cor(y,y_approx) > 0.9999

function func2(dd::Dict{Symbol,Float64})
    return sin(dd[:a])* cos(dd[:c])/ ((1+dd[:b])^2) +  sqrt(1.0+dd[:b])
end
func = func2
nodes  =  24
degree =  9
function_takes_Dict = true
limits = OrderedDict{Symbol,Tuple{Float64,Float64}}([:a, :b, :c] .=> [(-2.0,1.0), (0.1,0.15), (5.0,11.0)])
approxim = create_chebyshev_approximation(func, nodes, degree, limits, function_takes_Dict)
st =  0.6
d = vcat(Iterators.product(convert(Array{Float64,1}, limits[:a][1]:st:limits[:a][2]),convert(Array{Float64,1}, limits[:b][1]:st:limits[:b][2]),convert(Array{Float64,1}, limits[:c][1]:st:limits[:c][2]))...)
dd = DataFrame()
dd[:a], dd[:b], dd[:c]  = vcat.(d...)
y = Array{Float64,1}(undef, size(dd)[1])
for i in 1:length(y)
    dic = Dict{Symbol,Float64}([:a,:b,:c] .=> [dd[i,:a], dd[i,:b], dd[i,:c]])
    y[i] = func(dic)
end
y_approx = evaluate(approxim, dd)
y .- y_approx
maximum(abs.(y .- y_approx)) < 0.001
#using Distributions
#cor(y,y_approx) > 0.999999


text:
using DiDySim
using Base.Test

@testset "DiDySim.jl" begin
    # Write your own tests here.
    @test 1 == 2
end


text:
module GtkMarkdownTextView

    using Gtk
    import Gtk.GtkTextIter

    using Markdown

    export MarkdownTextView, MarkdownColors

    struct MarkdownColors
        font_size::Int
        color::String
        background::String
        highlight_color::String
        highlight_background::String
    end

    MarkdownColors() =  MarkdownColors(13, "#000", "#fff", "#111", "#eee")

    mutable struct MarkdownTextView <: GtkTextView

        handle::Ptr{Gtk.GObject}
        view::GtkTextView
        buffer::GtkTextBuffer

        function MarkdownTextView(m::Markdown.MD, prelude::String, mc::MarkdownColors = MarkdownColors())

            buffer = GtkTextBuffer()
            buffer.text[String] = prelude
            view = GtkTextView(buffer)

            style_css(view,
                "window, view, textview, buffer, text {
                    background-color: $(mc.background);
                    color: $(mc.color);
                    font-family: Monaco, Consolas, Courier, monospace;
                    margin:0px;
                }"
            )

            #set_gtk_property!(view, :margin_left, 1)
            view.monospace[Bool] = true
            view.wrap_mode[Bool] = true

            fs = mc.font_size

            Gtk.create_tag(buffer, "normal",    font = "$fs")
            Gtk.create_tag(buffer, "h1",        font = "bold $(fs+3)")
            Gtk.create_tag(buffer, "h2",        font = "bold $(fs+2)")
            Gtk.create_tag(buffer, "h3",        font = "bold $(fs+1)")
            Gtk.create_tag(buffer, "h4",        font = "bold $(fs)")
            Gtk.create_tag(buffer, "h5",        font = "$(fs)")
            Gtk.create_tag(buffer, "h6",        font = "$(fs-1)")
            Gtk.create_tag(buffer, "bold",      font = "bold $(fs)")
            Gtk.create_tag(buffer, "italic",    font = "italic $fs")
            Gtk.create_tag(buffer, "code",      font = "bold $fs",
                foreground=mc.highlight_color, background=mc.highlight_background)

            insert_MD!(buffer, m)
#            tag(buffer, "normal", 1, length(buffer))

            n = new(view.handle, view, buffer)
            Gtk.gobject_move_ref(n, view)
        end

        MarkdownTextView(m::String) = MarkdownTextView(Markdown.parse(m), "")
        MarkdownTextView(m::String, prelude::String, mc::MarkdownColors = MarkdownColors()) = MarkdownTextView(Markdown.parse(m), prelude, mc)
        MarkdownTextView(m::String, mc::MarkdownColors) = MarkdownTextView(Markdown.parse(m), "", mc)

    end

    function tag(buffer, what, i, j)
        Gtk.apply_tag(buffer, what,
            GtkTextIter(buffer, i), GtkTextIter(buffer, j)
        )
    end

    function style_css(w::Gtk.GtkWidget, css::String)
        sc = Gtk.G_.style_context(w)
        push!(sc, GtkCssProvider(data=css), 600)
    end

    function insert_MD!(buffer, m::Markdown.Header{N}, i) where N
        ip = i

        insert!(buffer, "    ")
        i += 4
        for el in m.text
            i = insert_MD!(buffer, el, i)
        end
        tag(buffer, "h$(min(N,4))", ip, i)
        i
    end

    function insert_MD!(buffer, m::Markdown.BlockQuote, i)
        insert!(buffer, "│  ")
        i += 3
        for el in m.content
            i = insert_MD!(buffer, el, i)
        end
        i
    end

    function insert_MD!(buffer, m::String, i)
        insert!(buffer, m)
        i += length(m)
    end

    function insert_MD!(buffer, m::Markdown.LaTeX, i)
        i = insert_MD!(buffer, m.formula, i)
    end

    function insert_MD!(buffer, m::Markdown.Paragraph, i)
    #    insert!(buffer, "\n\n")
    #    i += 2
        for el in m.content
            i = insert_MD!(buffer, el, i)
        end
        i
    end

    function insert_MD!(buffer, m::Markdown.Code, i)
        insert!(buffer, m.code)
        tag(buffer, "code", i, i+sizeof(m.code))
        i += length(m.code)
    end

    function insert_MD!(buffer, m::Markdown.List, i)

        marker = k -> m.ordered == -1 ? "•" : "$(k)."
        for (k, it) in enumerate(m.items)
            insert!(buffer, "    $(marker(k)) ")
            i += 6 + (m.ordered == 1)
            for el in it
                i = insert_MD!(buffer, el, i)
            end
            insert!(buffer, "\n")
            i += 1
        end
        i
    end

    tagname(m::Markdown.Italic) = "italic"
    tagname(m::Markdown.Bold) = "bold"

    function insert_MD!(buffer, m::T, i) where T <: Union{Markdown.Italic, Markdown.Bold}
        ip = i
        for el in m.text
            i = insert_MD!(buffer, el, i)
        end
        tag(buffer, tagname(m), ip, i)
        i
    end

    function insert_MD!(buffer, m, i)
        if isdefined(m, :text)
            for el in m.text
                i = insert_MD!(buffer, el, i)
            end
        end
        if isdefined(m, :content)
            for el in m.content
                i = insert_MD!(buffer, el, i)
            end
        end
        i
    end

    function insert_MD!(buffer, m::Markdown.MD)
        i = length(buffer)+1
        for el in m.content
            i = insert_MD!(buffer, el, i)
            insert!(buffer, "\n\n")
            i += 2
        end
    end

end


text:
export layout, subplot

"""
    layout(f, anchors)

Describes a custom grid layout allowing to specify where the axes origin should
be and how big (all relative to width/height of the figure). The `anchors` matrix
is a `n × 4` matrix where `n` is the total number of axes. The first two columns
indicate the relative position of the origin (so that the actual position is
`(anchors[k,1]W, anchors[k,2]H)` where `W,H=gcf().size`). The last two
columns indicate the size of the axes.

# Examples
```jl
julia> layout(f, [0.1 0.1 0.3 0.3;  # origin = (0.1W,0.1H)
                   0.1 0.5 0.3 0.3]) # size = (0.3W, 0.3H)
```
"""
function layout(f::Figure, anchors::Matrix{Float64})
    #  ______________________
    # |                      |
    # |  X    Y    W    H    |
    # |______________________|
    size(anchors, 2) == 4 || throw(ArgumentError("anchors must be of size ((nrows*ncols) × 4)"))
    all(0 .<= anchors .<= 1) || throw(ArgumentError("layout relative anchors must be between 0 and 1"))

    erase!(f)
    W, H = f.size
    # fill with Axes2D, if later there are axes3D it will just replace
    for i ∈ 1:size(anchors, 1)
        add_axes!(f, Axes2D{B}(parent=f.id,
                               origin=(anchors[i, 1]*W, anchors[i, 2]*H),
                               size=(anchors[i, 3]*W, anchors[i, 4]*H)))
    end
    return PreviewFigure(gcf())
end

"""
    subplot(a, b, c)
    subplot(abc)

Describe an automatic grid layout of axes of size `a × b` and selects axes `c`.
The axes are counted from the top left (1) increasing towards the right and the bottom.

# Examples
```jl
julia> subplot(222) # `2×2` grid layout, select top-right axes
julia> subplot(224) # selects bottom-right axes.
```

See also: [`layout`](@ref) to specify a custom layout.
"""
function subplot(nrows::Int, ncols::Int, idx::Int)
    1 <= nrows <= 9 || throw(ArgumentError("nrows must be between 1 and 9"))
    1 <= ncols <= 9 || throw(ArgumentError("ncols must be between 1 and 9"))
    1 <= idx <= nrows*ncols || throw(ArgumentError("idx must be between 1 and $(nrows*ncols)"))

    f = gcf()
    W, H = f.size
    wh_ratio = W/H
    # 1. if there are no axes in the current figure, add them
    if isempty(f.axes)
        # default layout is a grid
        #  ______________________
        # |                      |
        # | a1     a2     a3     |
        # |                      |
        # |                      |
        # | a4     a5     a6     |
        # |______________________|

        # this has been hand-tuned to give acceptable results for
        # (1, 2, 3) x (1, 2, 3). Beyond that is likely too cramped

        v_rel = 0.2 * wh_ratio * (1 + 0.2*(nrows-1))  # relative vertical margin
        vbox  = 1.0 / (nrows * (1.0 + v_rel) + v_rel) # relative vertical box size
        h_rel = 0.2 * (1 + 0.2*(ncols-1))             # relative horizontal margin
        hbox  = 1.0 / (ncols * (1.0 + h_rel) + h_rel) # relative horizontal boxsize

        vbox_m  = vbox * v_rel
        hbox_m  = hbox * h_rel
        vbox_wm = vbox + vbox_m
        hbox_wm = hbox + hbox_m
        hard_bottom = 0.01 * wh_ratio
        hard_left   = 0.01

        grid = zeros(nrows*ncols, 4)
        k = 1
        for r ∈ 1:nrows
            for c ∈ 1:ncols
                # increasing left to right
                # decreasing top to bottom
                grid[k, 1] = hard_left + hbox_m + (hbox_wm * (c-1)) # increasing left to right
                grid[k, 3] = hbox
                grid[k, 2] = hard_bottom + 1.0 - vbox_wm * r
                grid[k, 4] = vbox
                k += 1
            end
        end
        layout(f, grid)
    # 2. if there are axes, check that it matches, if it doesn't
    else
        length(f.axes) == nrows*ncols || throw(ArgumentError("the layout description does not " *
                              "match the current axes. If you want to change the layout of " *
                              "the current figure use `clf()` first to remove the existing axes."))
    end
    # 3. select the relevant axes and make them the current ones
    curax = f.axes[idx]
    GP_ENV["CUR_AXES"] = curax
    return nothing
end

function subplot(d::Int)
    111 <= d <= 999 || throw(ArgumentError("invalid description for the layout"))
    nrows = div(d, 100)
    dd    = d - 100nrows
    ncols = div(dd, 10)
    idx   = dd - 10ncols
    subplot(nrows, ncols, idx)
end


text:
#################### MCMC Simulation Engine ####################

"""
    mcmc(mc::ModelChains, iters::Integer; verbose::Bool=true, trees::Bool=false)

This function simulates additional draws from a model.

* `mc` is the results of a previous call to the mcmc function.

* `iters` indicates the number of draws to simulate.

* `verbose` controls whether to print progress statements to the console.

* `trees` indicates if the states of the model nodes describing tree structures should be stored as well.
"""
function mcmc(
    mc::ModelChains,
    iters::Integer;
    verbose::Bool = true,
    trees::Bool = false,
)::ModelChains

    thin = step(mc)
    last(mc) == div(mc.model.iter, thin) * thin ||
        throw(ArgumentError("chain is missing its last iteration"))

    mm = deepcopy(mc.model)
    mc2 = mcmc_master!(
        mm,
        mm.iter .+ (1:iters),
        mc.sim_params,
        burnin = last(mc),
        conv_storage = mc.conv_storage,
        samplers = mc.samplers,
    )

    if mc2.names != mc.names
        mc2 = mc2[:, mc.names, :]
    end
    ModelChains(
        vcat(mc, mc2),
        mc2.model,
        cat(mc.stats, mc2.stats, dims = 1),
        mc.stat_names,
        mc.sim_params,
        mc2.conv_storage,
        mc2.samplers,
    )
end

"""
    mcmc(m::Model, inputs::Dict{Symbol},
              inits::Vector{V} where V<:Dict{Symbol},
              iters::Integer; burnin::Integer=0, thin::Integer=1,
              chains::Integer=1, verbose::Bool=true, trees::Bool=false,
              params::SimulationParameters=SimulationParameters())

Simulate MCMC draws from the model `m`.

* `inputs` is a dictionary storing the values for input model nodes. Dictionary keys and values should be given for each input node.

* `inits` contains dictionaries with initial values for stochastic model nodes. Dictionary keys and values should be given for each stochastic node. Consecutive runs of the simulator will iterate through the vector’s dictionary elements.

* `iter` Specifies the number of draws to generate for each simulation run

* `burnin` specifies the number of initial draws to discard as a burn-in sequence to allow for convergence.

* `thin` is the step-size between draws to store in the output.

* `chains` specifies the number of simulation runs to perform.

* `verbose` indicates whether the sampler progress should be printed at the console.

* `trees` indicates if the states of the model nodes describing tree structures should be stored as well.

* `params` pass one Struct to set all simulation parameters, instead of setting each individually
"""
function mcmc(
    m::Model,
    inputs::Dict{Symbol},
    inits::Vector{V} where {V<:Dict{Symbol}},
    iters::Integer;
    burnin::Integer = 0,
    thin::Integer = 1,
    chains::Integer = 1,
    verbose::Bool = true,
    trees::Bool = false,
    params::SimulationParameters = SimulationParameters(),
)::ModelChains
    burnin = burnin == 0 ? params.burnin : burnin
    thin = thin == 1 ? params.thin : thin
    chains = chains == 1 ? params.chains : chains
    verbose = verbose == true ? params.verbose : verbose
    trees = trees == false ? params.trees : trees
    params = SimulationParameters(
        burnin,
        thin,
        chains,
        verbose,
        trees,
        params.asdsf,
        params.freq,
        params.min_splits,
    )

    params.asdsf &&
        !params.trees &&
        throw(ArgumentError("ASDSF can not be calculated without trees"))
    params.asdsf &&
        chains < 2 &&
        throw(ArgumentError("ASDSF can not be calculated one just one chain"))
    iters > params.burnin ||
        throw(ArgumentError("burnin is greater than or equal to iters"))
    length(inits) >= params.chains ||
        throw(ArgumentError("fewer initial values than chains"))

    mm::Model = deepcopy(m)
    setinputs!(mm, inputs)
    setinits!(mm, inits[1:params.chains])
    mm.burnin = burnin
    mcmc_master!(mm, 1:iters, params)
end

"""
  mcmc_master!(m::Model, window::UnitRange{Int}, sp::SimulationParameters;
               burnin::Int64=-1)::ModelChains

--- INTERNAL ---
Dispatchs parameters to corresponding functions to start the simulation and
builds the ModelChain object from the chains it receives.
"""
function mcmc_master!(
    m::Model,
    window::UnitRange{Int},
    sp::SimulationParameters;
    burnin::Int64 = -1,
    conv_storage::Union{Nothing,ConvergenceStorage} = nothing,
    samplers::Union{Nothing,Vector{Vector{Sampler}}} = nothing,
)::ModelChains

    chains = 1:sp.chains
    N = length(window)
    K = length(chains)

    sp.verbose &&
        println("MCMC Simulation of $N Iterations x $K Chain" * "s"^(K > 1) * "...")
    # this is necessary for additional draws from a model
    burnin = burnin == -1 ? sp.burnin : burnin

    states::Vector{ModelState} = m.states
    m.states = ModelState[]

    lsts =
        [Any[m, states[k], window, burnin, sp.thin, sp.trees, sp.verbose] for k in chains]
    if !isnothing(samplers)
        for k in chains
            lsts[k][1].samplers = samplers[k]
        end
    end

    results::Vector{Tuple{Chains,Model,ModelState}},
    stats::Array{Float64,2},
    statnames::Vector{AbstractString},
    conv_storage::Union{Nothing,ConvergenceStorage} =
        assign_mcmc_work(mcmc_or_convergence, lsts, sp, conv_storage)

    sims::Array{Chains} = Chains[results[k][1] for k = 1:K]
    model::Model = results[1][2]
    samplers = [res[2].samplers for res in results]

    model.states = ModelState[results[k][3] for k in sortperm(chains)]
    ModelChains(cat(sims..., dims = 3), model, stats, statnames, sp, conv_storage, samplers)
end

"""
  mcmc_worker!(args::Vector, ASDSF_step::Int64=0,
               rc::Union{Nothing, RemoteChannel}=nothing
               )::Tuple{Chains, Model, ModelState}

--- INTERNAL ---
Each call to this function computes a chain for a ModelChains Object.
"""
function mcmc_worker!(
    args::AbstractArray,
    ASDSF_step::Int64 = 0,
    rc::Union{Nothing,RemoteChannel} = nothing,
)::Tuple{Chains,Model,ModelState}
    m::Model,
    state::ModelState,
    window::UnitRange{Int},
    burnin::Integer,
    thin::Integer,
    store_trees::Bool,
    verbose::Bool,
    channel::Tuple{RemoteChannel,Int} = args
    llname::AbstractString = "likelihood"
    treeind::Int64 = 1
    m.iter = first(window) - 1

    relist!(m, state.value)
    initialize_samplers!(m)
    settune!(m, state.tune)
    pnames = vcat(names(m, true), llname)
    treenodes = Symbol[]
    for i in m.nodes
        if isa(i[2], Stochastic{<:GeneralNode})
            push!(treenodes, i[1])
        end
    end

    sim = Chains(
        last(window),
        length(pnames),
        start = burnin + thin,
        thin = thin,
        names = pnames,
        ntrees = length(treenodes),
        tree_names = treenodes,
    )

    for i in window

        sample!(m)
        if i > burnin
            if (i - burnin) % thin == 0
                sim[i, :, 1] = unlist(m, true)
                if store_trees
                    for (ind, tree_node) in enumerate(treenodes)
                        sim.trees[treeind, ind, 1] = newick(m[tree_node].value)
                    end # for
                    treeind += 1
                end # if
            end # if
            if !isnothing(rc) && (i - burnin) % ASDSF_step == 0
                trees::Vector{AbstractString} = []
                for (ind, tree_node) in enumerate(treenodes)
                    push!(trees, newick(m[tree_node].value))
                end # for
                # send trees to a RemoteChannel, so that convergence statistics can be calculated on a different worker
                put!(rc, trees)
            end # if
        end # if
        # send update to RemoteChannel --> while loop in logpdf function updates the ProgressMeter of this chain
        put!(channel[1], channel[2])
    end # for
    # signal to the assign_mcmc_work function, that this chain is finished
    put!(channel[1], -1)
    (sim, m, ModelState(unlist(m), gettune(m)))
end # mcmc_worker!

"""
  track(i::Integer, burnin::Integer, thin::Integer, sim::Chains, m::Model,
        store_trees::Bool, treeind::Integer, treenode)

--- INTERNAL ---
"""
function track(
    i::Integer,
    burnin::Integer,
    thin::Integer,
    sim::Chains,
    m::Model,
    store_trees::Bool,
    treeind::Integer,
    treenode,
)
    if i > burnin && (i - burnin) % thin == 0
        sim[i, :, 1] = unlist(m, true)
        if store_trees
            sim.trees[treeind, 1, 1] = newick(m[treenode].value)
            treeind += 1
        end # if
    end # if
end # track

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
"""
    module Simulate_FallingBall3

Models and simulates a sphere that is falling downwards as 6 dof body
(described by quaternions), **without** collision handling.

See also models

- [`Simulate_FallingBall1`](@ref) (falling ball constrained with prismatic joint),
- [`Simulate_FallingBall2`](@ref) (falling ball constrained with prismatic joint and collision handling).
"""
module Simulate_FallingBall3

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
ballMaterial   = Modia3D.Material(color="LightBlue", transparency=0.5)

@assembly FallingBall3(;h_start = 1.0, m = 0.1) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5))
   ground = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(1.0,0.1,1.0), nothing, groundMaterial),
                                                  r = [0.5, -0.05, 0.0])
   ball   = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(0.1), "Aluminium", ballMaterial),
                             fixed=false, r=[0.5,h_start,0.0])
end

fallingBall = FallingBall3()

# Modia3D.visualizeAssembly!( fallingBall )

simulationModel = Modia3D.SimulationModel( fallingBall, stopTime=0.6 )
result          = ModiaMath.simulate!(simulationModel, log=true)

ModiaMath.plot(result, ["ball.r[2]"])

println("... success of Simulate_FallingBall3.jl!")

end


text:
#=
quantum_bayesian:
- Julia version:
- Author: marcin
- Date: 2018-06-02
=#

const DiscreteQCPD = DiscreteMCPD{HermitianMatrix}

const AcausalStructure = BayesNet{DiscreteQCPD}

AcausalStructure() = BayesNet(DiscreteQCPD)


text:
using ITensors
using ITensorVisualizationBase

N = 10
sites(n) = Index([QN("Sz", 0) => 1, QN("Sz", 1) => 1]; tags="S=1/2,Site,n=$n")
l(n) = Index([QN("Sz", 0) => 10, QN("Sz", 1) => 10]; tags="Link,l=$n")
h(n) = Index([QN("Sz", 0) => 5, QN("Sz", 1) => 5]; tags="ham,Link,l=$n")
s⃗ = [sites(n) for n in 1:N]
l⃗ = [l(n) for n in 1:(N - 1)]
h⃗ = [h(n) for n in 1:(N - 1)]

# Add some more indices between two of the tensors
x = Index([QN("Sz", 0) => 2]; tags="X")
y = Index([QN("Sz", 0) => 2]; tags="Y")

n = 2
ψn1n2 = randomITensor(l⃗[n - 1], s⃗[n], s⃗[n + 1], l⃗[n + 1], dag(x), dag(y))
hn1 = randomITensor(dag(h⃗[n - 1]), s⃗[n]', dag(s⃗[n]), h⃗[n], x, y)
hn2 = randomITensor(dag(h⃗[n]), s⃗[n + 1]', dag(s⃗[n + 1]), h⃗[n + 1])
ELn0 = randomITensor(l⃗[n - 1]', h⃗[n - 1], dag(l⃗[n - 1]))
ERn2 = randomITensor(l⃗[n + 1]', dag(h⃗[n + 1]), dag(l⃗[n + 1]))

edge_labels = (; plevs=true)

R = @visualize fig1 ELn0 * ψn1n2 * hn1 * hn2 * ERn2 edge_labels = edge_labels vertex_size =
  50
@show R ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2

# Split it up into multiple contractions
R1 = @visualize fig2 ELn0 * ψn1n2 * hn1 edge_labels = edge_labels vertex_size = 50
R2 = @visualize fig3 R1 * hn2 * ERn2 edge_labels = edge_labels vertex_size = 50
@show R2 ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2

fig1, fig2, fig3


text:
"""
---
Apply SCFDMA demodulator to input signal and returns the T/F QAM matrix
# --- Syntax
qamRx	= scfdmaSigDecode(sigId,nFFT,nCP,allocSubcarriers,sizeDFT);
# --- Input parameters
# ---  Input parameters
- sigId	  : SCFDMA signal in time domain [Array{Complex{Float64},nbEch}]
- nFFT	  : FFT size [Int]
- nCp	  : Cyclic prefix size (in samples) [Int]
- allocSubcarriers : Vector of index of allocated subcarriers [Array{Int,nbSubcarriers}]
- sizeDFT	: Post processing size (DFT size)
- doPostProcessing : Do the post processing stage (IDFT): Default 1. In some case, we need the raw data (i.e data before postprocessing stage) such as in channel equalisation. In this case, the post processing should be done in a separate step (with the call of scfdmaPostProcessing function)
# ---  Output parameters
- qamRx  : Time frequency matrix : [Array{Complex{Float64},nbSubcarriers,nbSymb}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function scfdmaSigDecode(sigChan,nFFT,nCP,allocSubcarriers,sizeDFT,doPostProcessing=1)
	# Set 4 core for FFT computation
	#FFTW.set_num_threads(4)
	# --- Getting symbol size
	sizeSymb	   = nFFT + nCP;
	# --- Getting number of symbol to decode
	nbSymbDecode   = Int(floor(length(sigChan)/sizeSymb));
	# --- Initiate Rx T/F matrix
	qamRx		   = zeros(Complex{Float64},length(allocSubcarriers),nbSymbDecode);
	# --- SCFDMA demodulator
	for n = 1  : 1 : nbSymbDecode
		# --- Getting current symbol without CP
		currSymb  = sigChan[sizeSymb*(n-1).+collect(1+nCP:sizeSymb)];
		# --- TF
		fftEl	  = fft(currSymb,1);
		# --- Mapping in output matrix
		qamRx[:,n]		= fftEl[allocSubcarriers];
	end
	if Bool(doPostProcessing)
		qamRx	= scfdmaPostProcessing(qamRx,sizeDFT);
	end
	return qamRx;
end

"""
---
Apply Post processing for SC-FDMA (i.e DFT post-processing stage). This function can be called if scfdmaSigDecode is called with post processing flag to 0 (if frequency egalisation is done for example)
# --- Syntax
qamPost	  = scfdmaPostProcessing(qamMat,sizeDFT)
# --- Input parameters
- qamMat  : T/F matrix after Rx FFT (and before Rx IDFT) [Array{Complex{Float64}},nbSubcarriers,nbSymb]
- sizeDFT	: Precoder bloc size (often 12)
# --- Output parameters
- qamPost	: T/F matrix after IDFT processing [Array{Complex{Float64},nbSubcarriers,nbSymb}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function scfdmaPostProcessing(qamMat,sizeDFT)
	nbSymb		= size(qamMat,2);
	nbSubcarriers = size(qamMat,1);
	nbCarriers	= nbSubcarriers ÷ sizeDFT;
	qamIDFT		= zeros(Complex{Float64},nbSubcarriers,nbSymb);
	for iN = 1 : 1 : nbSymb
		for iC = 1 : 1 : nbCarriers
			# --- Apply DFT to current carrier, at current symbol
			qamIDFT[(iC-1)*sizeDFT .+ (1:sizeDFT),iN] = ifft(qamMat[(iC-1)*sizeDFT  .+ (1:sizeDFT),iN]);
		end
	end
	return qamIDFT
end

"""
---
Apply SCFDMA demodulator to input signal and returns the T/F QAM matrix
# --- Syntax
qamRx	= scfdmaSigDecode(sigId,nFFT,nCP,allocSubcarriers,sizeDFT);
# --- Input parameters
# ---  Input parameters
- sigId	  : SCFDMA signal in time domain [Array{Complex{Float64},nbEch}]
- scfdma  : SCFDMA structure [StrucSCFDMA]
- doPostProcessing : Do the post processing stage (IDFT): Default 1. In some case, we need the raw data (i.e data before postprocessing stage) such as in channel equalisation. In this case, the post processing should be done in a separate step (with the call of scfdmaPostProcessing function)
# ---  Output parameters
- qamRx  : Time frequency matrix : [Array{Complex{Float64},nbSubcarriers,nbSymb}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function scfdmaSigDecode(sigChan,scfdma::StrucSCFDMA,doPostProcessing=1)
	return scfdmaSigDecode(sigChan,scfdma.nFFT,scfdma.nCP,scfdma.allocatedSubcarriers,scfdma.sizeDFT,doPostProcessing);
end


text:
__precompile__(false)
module UniversalTensorBoard

using TensorFlow

Summaries = TensorFlow.summary

export set_tb_logdir, @tb_log, reset_tb_logs
export tensorboard

mutable struct TensorBoardSession
    logdir::String
    sess::Session
    writer::Summaries.FileWriter
    global_step::Int
end

function TensorBoardSession(logdir; overwrite=false)
    if overwrite
        rm(logdir; force=true, recursive=true)
    end
    mkpath(logdir)

    sess= Session(Graph())
    writer = Summaries.FileWriter(logdir; graph=sess.graph)
    TensorBoardSession(realpath(logdir), sess, writer, 0)
end

TensorBoardSession() = TensorBoardSession("tensorboard_logs", overwrite=true)
# normally the logs don't overwrite, but if you've not given a path, you clearly don't care.

const default_logging_session = Ref(TensorBoardSession())

####################################################################################################
# Interface

"""
    tensorboard()

Opens Tensorboard in the systems default browser.
"""
tensorboard(tb_sess=default_logging_session[]) = visualize(tb_sess.writer)

"""
    set_tb_logdir(logdir, overwrite=false)

Start a new log in the given directory
"""
function set_tb_logdir(logdir, overwrite=false)
    default_logging_session[] = TensorBoardSession(logdir, overwrite=overwrite)
end

"""
    reset_tb_logs()

Reset the current log, deleteing all information
"""
function reset_tb_logs()
    logdir = default_logging_session[].logdir
    default_logging_session[] = TensorBoardSession(logdir, overwrite=true)
end

macro tb_log(name)
    :(_tb_log($(esc(string(name))), $(esc(name))))
end

#################################################################
# Implmentation
function _tb_log(name, value; tb_sess = default_logging_session[])
    placeholder_node, summary_node  = get_nodes!(tb_sess, name, value)
    summaries = run(tb_sess.sess, summary_node,  Dict(placeholder_node=>value))
    write(tb_sess.writer, summaries, tb_sess.global_step)
    tb_sess.global_step+=1
    value
end

name_fix(name)=replace(name, r"[^A-Za-z0-9.]|\s", "_")
placeholder_node_name(name) = name_fix(name *"_input")
summary_node_name(name) = name_fix(name)

has_node(graph, name) = !isnull(TensorFlow.get_node_by_name(graph, name))

function get_nodes!(tb_sess, name, value)

    gg = tb_sess.sess.graph
    if has_node(gg, placeholder_node_name(name))
        (gg[placeholder_node_name(name)], gg[summary_node_name(name)])
    else
        TensorFlow.as_default(gg) do
            create_nodes!(name, value)
        end
    end
end

function create_nodes!(name, value::T) where T<:Number
    placeholder_node = placeholder(T; name=placeholder_node_name(name), shape=[])
    summary_node = Summaries.scalar(name, placeholder_node; name=summary_node_name(name))
    (placeholder_node, summary_node)
end

function create_nodes!(name, value::AbstractArray{T}) where T<:Number
    shape = collect(size(value))
    placeholder_node=placeholder(T; name=placeholder_node_name(name), shape=shape)
    summary_node = Summaries.histogram(name, placeholder_node; name=summary_node_name(name))
    (placeholder_node, summary_node)
end

end # modulee


text:
# A work in progress - beware, dragons!
# Feynman-Kleinert-1986-PRA
#
# Workbook to reimplement Feynman and Klenert's 1986 PRA "Effective classical partition function"

push!(LOAD_PATH,"../src/") # load module from local directory

using FeynmanKleinert

# https://doi.org/10.1103/PhysRevA.34.5080
#
# Effective classical partition functions. R. P. Feynman and H. Kleinert. Phys. Rev. A 34, 5080 â Published 1 December 1986
#
# Errata
#
# The form for the double well potential (page 34, RHS, third paragraph starting 'Another example is the double-well...'), should read
#
# $V(x)=-\frac{1}{2} x^2 + \frac{1}{4} g x^4 + \frac{1}{4g}$.
#
# Strangely this is correct in the captions of figure 2 and 3!

# Nb: Forked from Workbook by Jarv 2017-09-27

# Variables
# Potential energy g; either anharmonic strength or double-well setup
g=0.1976
# Thermodynamic Beta, natch
β=10

# Potential energy curves - V(x) functions
# Double well potential as given in figure 2 + 3 captions
# Hold my beer - I'm going to try metaprogramming
macro doublewell(g)
    return :(x -> -0.5*x^2 + 0.25*g*x^4 + 0.25/g)
end

# Anharmonic potential, following Kleinhert's book (p. 471); Also Figure 1 Feynman-Kleinert
macro anharmonicwell(g)
    return :(x^2/2+x^4*g/4)
end

# Harmonic oscillator - simple soln
macro harmonicwell(g)
    return :(0.5*x^2)
end

V=@doublewell(g)
println("@doublewell test. g=$g β=$β ")
println("V(0.0) = ",V(0.0))
println("Va2(0.0,1.0) = ",Va2(0.0,1.0,V))
println("Wtilde(0.0,4,4,V,β=$β) = ",Wtilde(0.0,1.0,4,V,β))

using Plots

# OK, let's have a look our actors:
# V, the bare potential
# Va2, the Gaussian-smoothed potential
# Wtilde, the Auxillary potential, including partition function magic / failure
#
xrange=-3:0.1:3
#
plot(x->V(x),xrange,label="V") # bare [potential]
#
for a2 in 0.5 #:0.1:0.5 # Gaussian smearing widths to apply
    plot!(x->Va2(x,a2,V)[1],xrange,label="Va2, a2=$a2")
    plot!(x->Wtilde(x,a2,5,V,β),xrange,label="Wtilde, a2=$a2") # bare [potential]
end

#plot!(ylim=(0,1)) # force display
plot!()

#### Print out a table of values for omega and A2 values ####
g=0.4
β=10

@printf("Wtilde(x0,a2,Ω,g=%f,β=%f)\n",g,β)

@printf("a2\\Ω")
for Ω in 0.1:0.2:2.0
    @printf(" %.3f",Ω)
end

for a2 in 0.1:0.1:2.0
    @printf("\n a2=%.3f",a2)
    for Ω in 0.1:0.2:2.0
        @printf(" %.3f",Wtilde(0.0,a2,Ω,V,β))
    end
end

# Replot fits from Feynman-Kleinert Table II
# - basically checking whether the Wtilde equation is reporting correctly,
# and whether (as described in the caption) only Wtilde for g=0.1976 has double well structure
@printf("\n\nFeynman-Kleinert Table II reads:\n g      E0=Wtilde\n 0.1976 0.650 \n 0.4    0.549 \n 4.0    0.598 \n 40.0   1.409\n")
β=3000
# Much higher than this, and results collapse to Inf. But essentially agrees with Table II now!

@printf("\nI calculate, β=%g: ",β)
# Wtilde(x0,a2,Ω)
xrange=-4:0.1:4
# These values, from Feynman and Kleinhert, table II
g=0.1976
@printf("\n%f %f",g,Wtilde(1.943,0.397,1.255,@doublewell(g),β))
plot(x->Wtilde(x,0.397,1.255,@doublewell(g),β),xrange,label="g=0.1976")

g=0.4
@printf("\n%f %f",g,Wtilde(0.0,1.030,0.486,@doublewell(g),β))
plot!(x->Wtilde(x,1.030,0.486,@doublewell(g),β),xrange,label="g=0.4")

g=4.0
@printf("\n%f %f",g,Wtilde(0.0,0.3059,1.634,@doublewell(g),β))
plot!(x->Wtilde(x,0.3059,1.634,@doublewell(g),β),xrange,label="g=4.0")

g=40.0
@printf("\n%f %f",g,Wtilde(0.0,0.1306,3.829,V,β))
plot!(x->Wtilde(x,0.1306,3.829,@doublewell(g),β),xrange,label="g=40.0")

plot!(ylim=(-2,10)) # force plot in notebook
@printf("\n")

### OK - let's go for the big one!
# Try and reproduce Fig 3; Double-well for varying g.

g=0.1976
β=1000

xrange=-2:0.05:2
#Ws=[ W(x,g,β) for x in xrange ]
#println(Ws)
#plot(Ws,xrange,label="W(x)")

plot(size=(400,500),fmt=:png) # force aspect ratio;
    # use PNG as changing size seems to break whatever the default is (SVG?)
for g in [.2, .3,.35,.4,.6] # range of g used in Fig 3. (look carefully at unbroken lines)
    plot!(x->W(x/sqrt(g),@doublewell(g),β), xrange, label="g=$g")
        # Note scaling of x-range within fn. argument
    plot!(x->@doublewell(g)(x/sqrt(g)), xrange, label="", linestyle = :dash)
        # Bare potential = infinite Temperature
end
plot!(ylim=(0,1.25)) # same as Fig 3 scale

savefig("figure3.png")


text:
## Multilink Attack algorithms

function projection(x::Float64,y::Float64,slope::Int)
  return y - slope*x
end

function mixedMLA(g::Network)
  breakingpoints, cuts = breakingPointsCuts(g)

  λ = breakingpoints[1][3]
  bounds = fill((typemax(Float64), 0.),λ) #(lower bound, upper bound) init

  k = λ
  for p in breakingpoints[2:end]
    if p[1] != -1.
      for i in (p[3]+1):(k-1)
        π = projection(p[1],p[2],i)
        A = min(p[2] - p[1]*p[3],bestMinkCut(g,cuts,i))
#        A = p[2] - p[1]*p[3]
        bounds[λ-i] = (π,A)
      end
      α = projection(p[1],p[2],p[3])
      bounds[λ-p[3]] = (α,α)
      k = p[3]
    end
  end
  return bounds
#  α0 = breakingpoints[end][2]
#  bounds[λ+1] = (α0,α0)
end

function successMLA(g::Network)
  bounds = mixedMLA(g)
  Σ = 0
  for i in bounds
    if i[1] == i[2]
      Σ += 1
    end
  end
  return Σ/length(bounds)
end


text:
using TomoForward
using XfromProjections

# img = imread("test_data/shepplogan512.png")[:,:,1]
# H, W = 128, 128
# img = imresize(img, H, W)

img = zeros(128, 128)
img[40:60, 40:60] .= 1.0
H, W = size(img)

nangles = 30
detcount = Int(floor(size(img,1)*1.4))
proj_geom = ProjGeom(1.0, detcount, LinRange(0,pi,nangles+1)[1:nangles])

# test line projection model
A = fp_op_parallel2d_line(proj_geom, size(img, 1), size(img, 2))
p = A * vec(img)
p = reshape(p, nangles, detcount)

niter=500
lambda = 0.01
u0 = zeros(size(img))
recon2d_tv_primaldual!(u0, A, p, niter, lambda)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
Base.@kwdef struct StructInformation
    typename_symbol::Symbol
    typename_string::String
    struct_supertype::String
    struct_field_names::Vector{String}
    struct_field_lines::Vector{String}
end

@inline function _to_string(struct_info::StructInformation)
    typename_symbol = struct_info.typename_symbol
    typename_string = struct_info.typename_string
    struct_supertype = struct_info.struct_supertype
    struct_field_names = struct_info.struct_field_lines
    struct_field_lines = struct_info.struct_field_lines

    lines = String[]
    push!(lines, "Base.@kwdef mutable struct $(typename_string) <: $(struct_supertype)")
    for struct_field_line in struct_field_lines
        push!(lines, "    $(struct_field_line)")
    end
    push!(lines, "end")
    push!(lines, "@inline StructTypes.StructType(::Type{T}) where {T <: $(typename_string)} = StructTypes.Mutable()")
    push!(lines, "# @inline StructTypes.omitempties(::Type{T}) where {T <: $(typename_string)} = true")
    result = strip(join(lines, "\n")) * "\n\n"
    return result
end