text:
using PowerModels, PowerPlots, Ipopt
using Colors:  @colorant_str

## Plot Network Status
mp = Dict{String,Any}(
    "active_line" => Dict(:color => :darkgreen, :size => 6),
    "inactive_line" => Dict(:color => :red, :size => 6),
    "active_bus" => Dict(:color => :green, :size => 15),
    "active_bus_load_shed" => Dict(:color => :orange, :size => 10),
    "inactive_bus" => Dict(:color => :red, :size => 15),
    "active_gen" => Dict(:color => :green, :size => 15),
    "inactive_gen" => Dict(:color => :red, :size => 15),
    "connector" => Dict(:color => colorant"lightgrey", :size => 2, :style => :dash)
)
case = PowerModels.parse_file("pglib_opf_case14_ieee.m")

inactive_devices = Dict("bus" =>  1:5, "gen" => 1:3, "branch" => 1:9)
for (item_type, ids) in inactive_devices
    for id in ids
        case[item_type]["$(id)"][pm_component_status[item_type]] = pm_component_status_inactive[item_type]
    end
end

plot_network_status(case; filename="network_status.png", plot_size=(800,800), membership_properties=mp)

## Plot Power Flow
mp=Dict{String,Any}(
            "max_power" => Dict(:color => colorant"red3", :size => 6),
            "min_power" => Dict(:color => colorant"black", :size => 6),
            "bus" => Dict(:color => colorant"green", :size => 12),
            "gen" => Dict(:color => colorant"blue", :size => 12),
            "connector" => Dict(:color => colorant"lightgrey", :size => 4, :style => :dash),
            "label" => Dict(:color => colorant"black", :size => 16, :fontfamily=>"Arial", :textalign=>:center, :offset => 0.1)
            )

case = PowerModels.parse_file("pglib_opf_case5_pjm.m")
result = run_opf(case, DCPPowerModel, Ipopt.Optimizer)
PowerModels.update_data!(case, result["solution"])

plot_power_flow(case; filename="power_flow.png", plot_size=(800,800), membership_properties=mp)

## Plot System Voltage
mp = Dict{String, Any}(
    "base_kv" => Dict(:palette => :Accent, :size=>[1,3]),
)

case = PowerModels.parse_file("pglib_opf_case118_ieee.m")
plot_system_voltage(case; filename="system_voltage.png", plot_size=(800,800), exclude_sources=true, membership_properties=mp)

## Pegase 1354
mp = Dict{String, Any}(
    "base_kv" => Dict(:palette => :Accent, :size=>[1,3]),
)

case = PowerModels.parse_file("pglib_opf_case1354_pegase.m")
plot_system_voltage(case; filename="pegase_1354.png", plot_size=(800,800), exclude_sources=true, membership_properties=mp)

## RTE 1888
mp = Dict{String, Any}(
    "base_kv" => Dict(:palette => :Accent, :size=>[1,3]),
)

case = PowerModels.parse_file("pglib_opf_case1888_rte.m")
plot_system_voltage(case; filename="rte_1888.png", plot_size=(800,800), exclude_sources=true, membership_properties=mp)


text:
using QMTK
using Compat.Test
using Compat.Iterators

@testset "Basic Interface" begin

    @test isperiodic(Chain{Fixed}(10)) == false
    @test isperiodic(Square{Fixed}(10, 10)) == false

    @test isperiodic(Chain{Periodic}(10)) == true
    @test isperiodic(Square{Periodic}(10, 10)) == true

    @test shape(Chain{Fixed}(10)) == (10, )
    @test shape(Square{Fixed}(10, 10)) == (10, 10)

    @test length(Chain{Fixed}(10)) == 10
    @test length(Square{Fixed}(10, 10)) == 100

end # testset

@testset "Chain Lattice" begin
    include("Chain.jl")
end

@testset "Square Lattice" begin
    include("Square.jl")
end

ext: jl


text:
#    Problem : GROUP A
#    *********
#    Another function with nontrivial groups and repetitious elements.
#
#    Original SIF Source:
#    N. Gould, private communication.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    COSINE.SIF classification OUR2-AN-V-0
#
#    Number of variables (variable)

f = (x) -> begin
    fx = 0.0
    for i in 1:lastindex(x)-1
        α = -0.5x[i+1]+x[i]^2
        fx += cos(α)
    end
    return fx
end

g! = (g, x) -> begin
    for i in 1:lastindex(x)-1
        α = -0.5x[i+1]+x[i]^2
        g[i] -= 2.0sin(α)x[i]
        g[i+1] += 0.5sin(α)
    end
    return g
end

fg! = (g, x) -> begin
    fx = 0.0
    for i in 1:lastindex(x)-1
        α = -0.5x[i+1]+x[i]^2
        fx += cos(α)
        g[i] -= 2.0sin(α)*x[i]
        g[i+1] += 0.5sin(α)
    end
    return fx, g
end

init = (n::Int=10000) -> begin
    n < 2 && @warn("COSINE: number of variables must be ≥ 2")
    return n, ones(n)
end

TestSet["COSINE"] = UncProgram("COSINE", f, g!, fg!, init)


text:
using Plots,ApproxFun;pyplot()

x = Fun(identity,[0.,10.])
f = sin(x^2)
g = cos(x)

h = f + g^2
r = roots(h)
rp = roots(differentiate(h))
plot(h;legend=false,grid=false)
scatter!(r,h(r))
scatter!(rp,h(rp))
xlabel!("\$x\$")
ylabel!("\$h(x)\$")
xlims!(0,10)
PyPlot.savefig("extrema.png",dpi=300)
println("First image done")

x = Fun(identity,[-1000.,200.])
d = domain(x)
D = Derivative(d)
B = dirichlet(d)
L = D^2 - x
u = [B;L] \ [airyai(d.a),airyai(d.b)]
plot(u;legend=false,grid=false)						    # Requires Gadfly or PyPlot
xlabel!("\$x\$");ylabel!("\${\\rm Ai}(x)\$")
PyPlot.savefig("airy.png",dpi=300)

println("Second image done")

d = Interval(-π,π)
a = Fun(t-> 1+sin(cos(2t)),d)
D = Derivative(d)
L = D + a
f = Fun(t->exp(sin(10t)),d)
B = periodic(d,0)
uChebyshev = [B;L]\[0.,f]

d = PeriodicSegment([-π,π])
a = Fun(t-> 1+sin(cos(2t)),d)
D = Derivative(d)
L = D + a
f = Fun(t->exp(sin(10t)),d)
uFourier = L\f

length(uFourier)/length(uChebyshev),2/π
plot(real(uFourier);legend=false,grid=false)						    # Requires Gadfly or PyPlot
xlims!(d.a,d.b);xlabel!("\$t\$");ylabel!("\$u(t)\$")
PyPlot.savefig("periodic.png",dpi=300)

println("Third image done")

f = abs(Fun(sin,[-5.,5.]))
d = domain(f)
x = ApproxFun.sample(f,10000)
plot(f/sum(f);legend=false,grid=false)                           # Requires Gadfly or PyPlot
histogram!(x;normed=true,nbins=100)
xlims!(leftendpoint(d),rightendpoint(d));ylims!(0.0,0.18);xlabel!("\$x\$");ylabel!("Density")
PyPlot.savefig("Sample.png",dpi=300)

println("Fourth image done")

## Nonlinear BVP

x=Fun()
u0=0x

N=u->[u(-1.)-1.,u(1.)+0.5,0.001u''+6*(1-x^2)*u'+u^2-1.]
u=newton(N,u0)
plot(u;legend=false,grid=false)
PyPlot.savefig("nbvp.png",dpi=300)

println("Fifth image done")

## Multivariate

d = ChebyshevInterval()^2                            # Defines a rectangle

u = [dirichlet(d);lap(d)+100I]\[[1.0,1.0,1.0,1.0]]      # First four entries of rhs are
                                            # boundary conditions
surface(u,grid=false)                                     # contour plot

PyPlot.savefig("helmholtz.png",dpi=300)


text:
rootdir = Pkg.dir("Dyno")
filename = joinpath(rootdir,"examples","rbc.mod")

import Dyno
@time model = Dyno.import_model(filename, print_code=false)

model.functions.f_static

y = [model.calibration[p] for p in model.symbols[:endogenous]]
e = [model.calibration[p] for p in model.symbols[:exogenous]]
p = [model.calibration[p] for p in model.symbols[:parameters]]
using Dolang
model.functions.f_static(Der{1},[y;e],p)
model.functions.f_dynamic(Der{1},[y;y;y;e],p)

exo = [0.0, 0.1, 0.05, 0.0]''

sim = Dyno.deterministic(model, exo)

dr = Dyno.solve(model)


text:
#!/usr/bin/env julia6
using MyStats.condmean
using ReadGlobal: readfield, readfield!
using JLD2

function main(condname::AbstractString, inpname::AbstractVector{<:AbstractString}, outfilename::AbstractVector{<:AbstractString})
    @load condname ind bins
    field = readfield(inpname[1])
    for i in linearindices(inpname)
        readfield!(inpname[i],field)
        r,err = condmean(field,ind)
        writedlm(outfilename[i],hcat(bins,r,err))
    end
    return nothing
end

function main(args::Vector{<:AbstractString})
    condname = args[1]
    iseven(length(args)) && error()
    p = div(length(args) - 1, 2) + 1
    inpname = args[2:p]
    outfilename = args[p+1:end]
    main(condname,inpname,outfilename)
end

main(ARGS)


text:
using Test
using PandaModels; const _PdM = PandaModels
import PowerModels; const _PM = PowerModels

_PM.silence()

pdm_path = joinpath(dirname(pathof(PandaModels)), "..")
data_path = joinpath(pdm_path, "test", "data")

case_pm =  joinpath(data_path, "test_pm.json")
case_pf = joinpath(data_path, "test_pf.json")
case_opf = joinpath(data_path, "test_opf.json")
case_custom = joinpath(data_path, "test_custom.json")
case_ots = joinpath(data_path, "test_ots.json")
case_tnep = joinpath(data_path, "test_tnep.json")

@testset "PandaModels.jl" begin
        @testset "test internal functions" begin
                pm = _PdM.load_pm_from_json(case_pm)

                @test length(pm["bus"]) == 6
                @test length(pm["gen"]) == 3
                @test length(pm["branch"]) == 7
                @test length(pm["load"]) == 3

                model =_PdM.get_model(pm["pm_model"])
                @test string(model) == "PowerModels.DCPPowerModel"

                solver = _PdM.get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
                pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

                @test string(solver.optimizer_constructor) == "Ipopt.Optimizer"

        end

        @testset "test exported executive functions" begin
                @testset "test for run_powermodels" begin
                        result=run_powermodels_opf(case_opf)
                        @test isa(result, Dict{String,Any})
                        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
                        @test isapprox(result["objective"], 144.85; atol = 1e0)
                        @test result["solve_time"] > 0.0
                end
                @testset "test for run_powermodels_powerflow" begin
                        result=run_powermodels_powerflow(case_pf)
                        @test isa(result, Dict{String,Any})
                        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
                        @test isapprox(result["objective"], 0.0; atol = 1e0)
                        @test result["solve_time"] > 0.0
                end
                @testset "test for powermodels_custom" begin
                        result=run_powermodels_custom(case_custom)
                        @test isa(result, Dict{String,Any})
                        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
                        @test isapprox(result["objective"], 144.85; atol = 1e0)
                        @test result["solve_time"] > 0.0
                end
                @testset "test for powermodels_tnep" begin
                        result=run_powermodels_tnep(case_tnep)
                        @test isa(result, Dict{String,Any})
                        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
                        @test isapprox(result["objective"], 0.0; atol = 1e0)
                        @test result["solve_time"] > 0.0
                end

                @testset "test for powermodels_ots" begin
                        result=run_powermodels_ots(case_ots)
                        @test isa(result, Dict{String,Any})
                        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
                        @test isapprox(result["objective"], 14810.0; atol = 100.0)
                        @test result["solve_time"] > 0.0
                end
        end
end


text:
"""
loads and do the main processing of data in arrays of intrest (padding of shmem will be analyzed separately)

"""
module ProcessMainData
using  StaticArrays,..CUDAGpuUtils ,..HFUtils, CUDA, ..ProcessPadding
export executeDataIterFirstPass,executeDataIterOtherPasses,processMaskData,executeDataIterFirstPassWithPadding

"""
   loads main values from analyzed array into shared memory and to locArr - which live in registers
"""

macro loadMainValues()
        @iter3dWithVal  dataBlockDims loopX loopY loopZ blockBeginingX blockBeginingY blockBeginingZ analyzedArr begin
        #val is given by macro as value of this x,y,z
        locArr|= val << (zIter-1)
        processMaskData( val, zIter, resShmem)
        #zIter given in macro as we are iterating in this spot
        sourceShmem[threadIdxX(), threadIdxY(), zIter]
    end
end #loadMainValues

"""
 validates data is of our intrest
"""
macro validateData()
    @iter3dW  dataBlockDims loopX loopY loopZ blockBeginingX blockBeginingY blockBeginingZ resShemVal begin
        locVal::Bool = @inbounds  (locArr>>(zIter-1) & 1)
        resShemVal::Bool = @inbounds resShmem[threadIdxX()+1,threadIdxY()+1,zIter+1]
        locValOrShmem = (locVal | resShemVal)
        #those needed to establish weather data block will remain active
        isMaskFull= locValOrShmem & isMaskFull
        @ifverr zzz isMaskEmpty = ~locValOrShmem & isMaskEmpty
        if(!locVal && resShemVal)
              innerValidate(analyzedArr,referenceArray,x,y,z,privateResArray,privateResCounter,iterationnumber,sourceShmem  )
        end#if
     end#3d iter

 end  #validateData

"""
this will be invoked when we know that we have a true in a spot that was false before this dilatation step and its task is to set to true appropriate spot in global array
- so proper dilatation
check weather we have true also in reference array - if so we  need to add this spot to the block result list in case we are invoke it from padding we need to look even futher into the
next block data to establish could this spot be activated from there
"""
  function innerValidate(analyzedArr,referenceArray,x,y,z,privateResArray,privateResCounter,iterationnumber,sourceShmem  )
            # setting value in global memory
            @inbounds  analyzedArr[x,y,z]= true
            # if we are here we have some voxel that was false in a primary mask and is becoming now true - if it is additionaly true in reference we need to add it to result

            if(@inbounds referenceArray[x,y,z])
                #results now are stored in a matrix where first 3 entries are x,y,z coordinates entry 4 is in which iteration we covered it and entry 5 from which direction - this will be used if needed
                #privateResCounter privateResArray are holding in metadata blocks results and counter how many results were already added
                #in each thread block we will have separate rescounter, and res array for goldboolpass and other pass
               direction=  @ifverr zzz  getDir(sourceShmem) | 0
               @append  privateResArray privateResCounter  [x,y,z,iterationnumber, direction]

            end#if
  end#innerValidate

"""
Help to establish should we validate the voxel - so if ok add to result set, update the main array etc
  in case we have some true in padding
  generally we need just to get idea if
    we already had true in this very spot - if so we ignore it
    can this spot be reached by other voxels from the block we are reaching into - in other words padding is analyzing the same data as other block is analyzing in its main part
      hence if the block that is doing it in main part will reach this spot on its own we will ignore value from padding

  in order to reduce sears direction by 1 it would be also beneficial to know from where we had came - from what direction the block we are spilled into padding
"""
function isPaddingValToBeValidated(dir,analyzedArr, x,y,z )::Bool

if(dir!=5)  if( @inbounds resShmem[threadIdxX(),threadIdxY(),zIter-1]) return false  end end #up
if(dir!=6)  if( @inbounds  resShmem[threadIdxX(),threadIdxY(),zIter+1]) return false  end  end #down

if(dir!=1)   if( @inbounds  resShmem[threadIdxX()-1,threadIdxY(),zIter]) return false  end  end #left
if(dir!=2)   if( @inbounds   resShmem[threadIdxX()+1,threadIdxY(),zIter]) return false  end end  #right

if(dir!=4)   if(  @inbounds  resShmem[threadIdxX(),threadIdxY()+1,zIter]) return false  end  end #front
if(dir!=3)   if( @inbounds  resShmem[threadIdxX(),threadIdxY()-1,zIter]) return false  end end  #back
  #will return true only in case there is nothing around
  return true
end

    """
    now in case we  want later to establish source of the data - would like to find the true distances  not taking the assumption of isometric voxels
    we need to store now data from what direction given voxel was activated what will later gratly simplify the task of finding the true distance
    we will record first found true voxel from each of six directions
                     top 6
                    bottom 5
                    left 2
                    right 1
                    anterior 3
                    posterior 4
    """
   function getDir(sourceShmem)
    if( @inbounds sourceShmem[threadIdxX(),threadIdxY(),zIter-1]) return 6  end  #up
    if( @inbounds  sourceShmem[threadIdxX(),threadIdxY(),zIter+1]) return 5  end #down

    if( @inbounds  sourceShmem[threadIdxX()-1,threadIdxY(),zIter]) return 2  end #left
    if( @inbounds   sourceShmem[threadIdxX()+1,threadIdxY(),zIter]) return 1  end #right

    if(  @inbounds  sourceShmem[threadIdxX(),threadIdxY()+1,zIter]) return 3  end #front
     if( @inbounds  sourceShmem[threadIdxX(),threadIdxY()-1,zIter]) return 4  end #back
  end#getDir

####to 3d iter data
#       ---- to iteration3d we can get a function from macro to a generalized function and then produce multiple macros that will be chosen  based on multiple dispatch
#          same will get the macro iter3d with val so value will also be available  also we need to   we can also specify bounds safe and not safe variant - what will reduce number of ifs
#                     - also we can consider
#     ----------- need to add in 3 dim iter a possibility to customize the way how we define offsets in x,y,z so instead of grid or block dim we need to have ability to do it separately as extra arguments
# ifverr - will return only this expression that is compatible with version number supplied
###

function executeDataIterWithPadding(analyzedArr, referenceArray,blockBeginingX
                                ,blockBeginingY,blockBeginingZ,resShmem,sourceShmem,resArray,resArraysCounter
                                ,currBlockX,currBlockY,currBlockZ,isPassGold,metaData,metadataDims
                                ,mainQuesCounter,mainWorkQueue,iterationNumber,debugArr, loopX,loopY,loopZ, dataBlockDims,
                                privateResCounter, blockResCounter)

    #some data cleaning
    locArr::UInt32 = UInt32(0)
    # locFloat::Float32 = Float32(0.0)
    isMaskFull::Bool= true
    isMaskOkForProcessing::Bool = true
    offset = 1
    ############## upload data
    @loadMainValues
    syncthreads()
    ########## check data aprat from padding
    #can be skipped if we have the block with already all results analyzed
    if(getIsTotalFPorFNnotYetCovered(resshmem )   )
        @validateData
    end
    #processing padding
    @processPadding

        if(threadIdxY()==5 && threadIdxX()==5 && (resShmem[2,2,6] || resShmem[2,2,7]))
            metaData[currBlockX,currBlockY,currBlockZ,isPassGold+1]=false # we set is inactive
        end#if
        if(threadIdxY()==6 && threadIdxX()==6 && (resShmem[2,2,6] || resShmem[2,2,7]))
            metaData[currBlockX,currBlockY,currBlockZ,isPassGold+3]=true # we set is as full
        end#if

end#executeDataIter

"""
we are processing padding
"""
macro processPadding()
    #so here we utilize iter3 with 1 dim fixed
    @unroll for  dim in 1:3, numb in [1,34]
        @iter3dFixed dim numb if( isPaddingValToBeValidated(dir,analyzedArr, x,y,z ))
            if(isPaddingToBeAnalyzed(resShmem,dim,numb ))
                innerValidate(analyzedArr,referenceArray,x,y,z,privateResArray,privateResCounter,iterationnumber,sourceShmem  )
                #   --- here we need also to set appropriate value in metadata marking that block in given direction marked as to be acivated from padding     all still need to check is th ere any block at all
                #              so check metadata dims
            end

        end#iter3dFixed
    end#for

end

# """
# specializes executeDataIterFirstPass as it do not consider possibility of block being empty
# """
# function executeDataIterOtherPasses(analyzedArr, refAray,iterationNumber ,blockBeginingX,blockBeginingY,blockBeginingZ,isMaskFull,resShmem,locArr,resArray,resArraysCounter)
#     @unroll for zIter in UInt16(1):32# most outer loop is responsible for z dimension
#         local locBool = testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]
#         locArr|= locBool << zIter
#        processMaskData( locBool, zIter, resShmem)
#        CUDA.unsafe_free!(locBool)
#     end#for
#     sync_threads() #we should have in resShmem what we need
#     @unroll for zIter in UInt16(1):32 # most outer loop is responsible for z dimension - importnant in this loop we ignore padding we will deal with it separately
#          local locBoolRegister = (locArr>>zIter & UInt32(1))==1
#          local locBoolShmem = resShmem[threadIdxX()+1,threadIdxY()+1,zIter+1]
#          validataDataOtherPass(locBoolRegister,locBoolShmem,isMaskFull,isMaskEmpty,blockBeginingX,blockBeginingY,blockBeginingZ,analyzedArr, refAray,resArray,resArraysCounter)
#        CUDA.unsafe_free!(locBool)

#     end#for
# end#executeDataIter

#numb>>1 & UInt32(1)

"""
uploaded data from shared memory in amask of intrest gets processed in this function so we need to
    - save it to registers (to locArr)
    - save to the 6 surrounding voxels in shared memory intermediate results
            - as we also have padding we generally start from spot 2,2 as up and to the left we have 1 padding
            - also we need to make sure that in corner cases we are getting to correct spot
"""
function processMaskData(maskBool::Bool
                         ,zIter::UInt8
                         ,resShmem
                          ) #::CUDA.CuRefValue{Int32}
    # save it to registers - we will need it later
    #locArr[zIter]=maskBool
    #now we are saving results evrywhere we are intrested in so around without diagonals (we use supremum norm instead of euclidean)
    #locArr.x|= maskBool << zIter
    if(maskBool)
        @inbounds resShmem[threadIdxX()+1,threadIdxY()+1,zIter]=true #up
        @inbounds  resShmem[threadIdxX()+1,threadIdxY()+1,zIter+2]=true #down

        @inbounds  resShmem[threadIdxX(),threadIdxY()+1,zIter+1]=true #left
        @inbounds   resShmem[threadIdxX()+2,threadIdxY()+1,zIter+1]=true #right

        @inbounds  resShmem[threadIdxX()+1,threadIdxY()+2,zIter+1]=true #front
        @inbounds  resShmem[threadIdxX()+1,threadIdxY(),zIter+1]=true #back
    end#if

end#processMaskData

"""
-so we uploaded all data that we consider new - around voxels that are "true"  but we can be sure that some of those were already true earlier
    possibly it can be marked also by some other neighbouring thread in this particular sweep
    in order to reduce writes to global memory we need to check with registers wheather it is alrerady in a mask - and we will write it to global memory only if it was not
    if the true is in shmem but not in register we write it to global memory - if futher it is also present in other mask (that we are comparing with now)
    we write it also to global result array
- updata isMaskFull and isMaskEmpty if needed using data from registers and shmem - so later we will know is this mask s full or empty
- we need to take special care for padding - and in case we would find anything there we need to mark appropriate neighbouring block to get activated
    save result if it did not occured in other mask and write it to global memory

locVal - value from registers
shmemVal - value associated with this thread from shared memory - where we marked neighbours ...
resShmem - shared memory with our preliminary results
isMaskFull, isMaskEmpty - register values needed to specify weather we have full or empty or neither block
blockBeginingX,blockBeginingY,blockBeginingZ - coordinates where our block is begining - will be used as offset by our threads
masktoUpdate - mask that we analyzed and now we write to data about dilatation
maskToCompare - the other mask that we need to check before we write to result array
resArray
        """
function validataDataFirstPass(locVal::Bool
                    ,shmemVal::Bool
                    ,isMaskFull#::CUDA.CuRefValue{Bool}
                    ,isMaskEmpty#::CUDA.CuRefValue{Bool}
                    ,blockBeginingX::UInt8
                    ,blockBeginingY::UInt8
                    ,blockBeginingZ::UInt8
                    ,maskToCompare
                    ,masktoUpdate
                    ,resArray
                    ,resArraysCounter
                    ,zIter::UInt8)::Bool
    #when this one and previous is true it will still be true

return

end

function IterToValidate()::Bool
    @unroll for zIteB::UInt8 in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
        # locBoolRegister::Bool = getLocalBoolRegister(locArr,zIter)
            #  locBoolShmem::Bool = getLocalBoolShemem(resShmem, zIter)
            # ProcessMainData.validataDataFirstPass(locBoolRegister,locBoolShmem,resShmem,isMaskFull,isMaskEmpty,blockBeginingX,blockBeginingY,blockBeginingZ,testArrInn, referenceArray,resArray,resArraysCounter,zIter)
            #CUDA.unsafe_free!(locBoolRegister)
            # CUDA.unsafe_free!(locBoolShmem)
         end#for
     return true
end

"""
specializes validataDataFirstPass ignoring case of potentially empty mask
iterationNumber - in which iteration we are currently - the bigger it is the higher housedorfrf,,

"""
function validataDataOtherPass(locVal::Bool
                                ,shmemVal::Bool
                                ,isMaskEmpty::MVector{1, Bool}
                                ,blockBeginingX,blockBeginingY,blockBeginingZ
                                ,maskToCompare
                                ,masktoUpdate
                                ,resArray
                                ,iterationNumber::UInt16
                                ,resArraysCounter
                                ,zIter)
    #when this one and previous is true it will still be true
    setIsFull!((locVal | shmemVal),isMaskEmpty  )
    if(!locVal && shmemVal)
    # setting value in global memory
    masktoUpdate[x,y,z+32]= true
    # if we are here we have some voxel that was false in a primary mask and is becoming now true - if it is additionaly true in reference we need to add it to result
    if(maskToCompare[x,y,z+32])
    resArray[x,y,z+32]=iterationNumber
    CUDA.atomic_inc!(pointer(resArraysCounter), UInt16(1))

    end#if
    end#if

end

end#ProcessMainData

# function setIsFull!(locValOrShmem::Bool
#     ,isMaskFull::MVector{1, Bool})
# isMaskFull[1]= locValOrShmem & isMaskFull[1]
# end#setIsFullOrEmpty

# """
# loads and do the main processing of data in arrays of intrest (padding of shmem will be analyzed separately)
# analyzedArr - array we are currently dilatating
# refAray -array we are referencing (we do not dilatate it only check against it )
# iterationNumber - at what iteration of dilatation we are - so how many dilatations we already performed
# blockBeginingX,blockBeginingY,blockBeginingZ - coordinates where our block is begining - will be used as offset by our threads
# isMaskFull,isMaskEmpty - enables later checking is mask is empty full or neither
# resShmem - shared memory bit array of the same dimensions as data block but plus 2
# sourceShmem - bit array of the same dimensions as data block
# locArr - local bit array of thread
# resArray- 3 dimensional array where we put results
# loopX,loopY,loopZ - how many times main loops need to be invoked in order for the thread block to cover all of  data block
# dataBlockDims - size of data block - so max x y and z possible index in those dimensions
# """
# function executeDataIterFirstPass(analyzedArr, referenceArray,blockBeginingX,blockBeginingY,blockBeginingZ,resShmem,sourceShmem,resArray,resArraysCounter, loopX,loopY,loopZ, dataBlockDims)
#     locArr = Int32(0)
#     isMaskFull= true
#     isMaskEmpty = true
#     #locArr.x |= true << UInt8(2)

#     ############## upload data
#     ---- to iteration3d we can get a function from macro to a generalized function and then produce multiple macros that will be chosen  based on multiple dispatch
#          same will get the macro iter3d with val so value will also be available  also we need to   we can also specify bounds safe and not safe variant - what will reduce number of ifs
#                     - also we can consider
#     ###step 1
#     @iter3dWithVal  dataBlockDims loopX loopY loopZ blockBeginingX blockBeginingY blockBeginingZ analyzedArr begin
#         #val is given by macro as value of this x,y,z
#         locArr|= val << (zIter-1)
#         processMaskData( val, zIter, resShmem)
#         #zIter given in macro as we are iterating in this spot
#         sourceShmem[threadIdxX(), threadIdxY(), zIter]
#     end
#     syncthreads()

#     ##step 2
#     ########## check data aprat from padding
#     @iter3dWithVal  dataBlockDims loopX loopY loopZ blockBeginingX blockBeginingY blockBeginingZ analyzedArr begin
#         locVal::Bool = @inbounds  (locArr>>(zIter-1) & 1)
#         resShemVal::Bool = @inbounds resShmem[threadIdxX()+1,threadIdxY()+1,zIter+1]
#     end

#     ########## check data aprat from padding
#     @unroll for zIter::UInt8 in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
#         locVal::Bool = @inbounds  (locArr>>(zIter-1) & 1)
#         shmemVal::Bool = @inbounds resShmem[threadIdxX()+1,threadIdxY()+1,zIter+1]
#         # CUDA.@cuprint "locVal $(locVal)  shmemVal $(shmemVal) \n  "
#         locValOrShmem = (locVal | shmemVal)
#         isMaskFull= locValOrShmem & isMaskFull
#         isMaskEmpty = ~locValOrShmem & isMaskEmpty

#         #CUDA.@cuprint "locVal $(locVal)  shmemVal $(shmemVal) \n  "
#         if(!locVal && shmemVal)
#             # setting value in global memory
#             @inbounds  analyzedArr[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]= true
#             # if we are here we have some voxel that was false in a primary mask and is becoming now true - if it is additionaly true in reference we need to add it to result
#             isInReferencaArr::Bool= @inbounds referenceArray[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]
#             if(isInReferencaArr)
#                 #CUDA.@cuprint "isInReferencaArr $(isInReferencaArr) \n  "
#                 @inbounds  resArray[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]=UInt16(1)
#                 #CUDA.atomic_inc!(pointer(resArraysCounter), Int32(1))
#                 atomicallyAddOneInt(resArraysCounter)
#             end#if
#         end#if
#       end#for

# end#executeDataIter


text:
using Test

let x, y, x
    x = TypeVar(:a)
    y = TypeVar(:a)
    z = TypeVar(:a)
    @test repr(UnionAll(z, UnionAll(x, UnionAll(y, Tuple{x,y,z})))) == "Tuple{a1, a2, a} where {a, a1, a2}"
    @test repr(UnionAll(z, UnionAll(x, UnionAll(y, Tuple{z,y,x})))) == "Tuple{a, a2, a1} where {a, a1, a2}"
end

let x = TypeVar(:_, Number), y = TypeVar(:_, Number)
    @test repr(UnionAll(x, UnionAll(y, Pair{x,y}))) == "Pair{_1, _2} where {_1<:Number, _2<:Number}"
    @test repr(UnionAll(y, UnionAll(x, Pair{x,y}))) == "Pair{_2, _1} where {_1<:Number, _2<:Number}"
    @test repr(UnionAll(x, UnionAll(y, Pair{UnionAll(x, Ref{x}),y}))) == "Pair{Ref{_1} where _1<:Number, _1} where _1<:Number"
    @test repr(UnionAll(y, UnionAll(x, Pair{UnionAll(y, Ref{x}),y}))) == "Pair{Ref{_2}, _1} where {_1<:Number, _2<:Number}"
end

@testset "matrix printing" begin
    # print_matrix should be able to handle small and large objects easily, test by
    # calling show. This also indirectly tests print_matrix_row, which
    # is used repeatedly by print_matrix.
    # This fits on screen:
    @test replstr(Matrix(1.0I, 10, 10)) == "10×10 Matrix{Float64}:\n 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0\n 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0\n 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0\n 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0\n 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0\n 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0\n 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0\n 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0\n 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0\n 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0"
    # an array too long vertically to fit on screen, and too long horizontally:
    @test replstr(Vector(1.0:100.0)) == "100-element Vector{Float64}:\n   1.0\n   2.0\n   3.0\n   4.0\n   5.0\n   6.0\n   7.0\n   8.0\n   9.0\n  10.0\n   ⋮\n  92.0\n  93.0\n  94.0\n  95.0\n  96.0\n  97.0\n  98.0\n  99.0\n 100.0"
    @test occursin(r"1×100 adjoint\(::Vector{Float64}\) with eltype Float64:\n 1.0  2.0  3.0  4.0  5.0  6.0  7.0  …  95.0  96.0  97.0  98.0  99.0  100.0", replstr(Vector(1.0:100.0)'))
    # too big in both directions to fit on screen:
    @test replstr((1.0:100.0) * (1:100)') == "100×100 Matrix{Float64}:\n   1.0    2.0    3.0    4.0    5.0    6.0  …    97.0    98.0    99.0    100.0\n   2.0    4.0    6.0    8.0   10.0   12.0      194.0   196.0   198.0    200.0\n   3.0    6.0    9.0   12.0   15.0   18.0      291.0   294.0   297.0    300.0\n   4.0    8.0   12.0   16.0   20.0   24.0      388.0   392.0   396.0    400.0\n   5.0   10.0   15.0   20.0   25.0   30.0      485.0   490.0   495.0    500.0\n   6.0   12.0   18.0   24.0   30.0   36.0  …   582.0   588.0   594.0    600.0\n   7.0   14.0   21.0   28.0   35.0   42.0      679.0   686.0   693.0    700.0\n   8.0   16.0   24.0   32.0   40.0   48.0      776.0   784.0   792.0    800.0\n   9.0   18.0   27.0   36.0   45.0   54.0      873.0   882.0   891.0    900.0\n  10.0   20.0   30.0   40.0   50.0   60.0      970.0   980.0   990.0   1000.0\n   ⋮                                  ⋮    ⋱                          \n  92.0  184.0  276.0  368.0  460.0  552.0     8924.0  9016.0  9108.0   9200.0\n  93.0  186.0  279.0  372.0  465.0  558.0     9021.0  9114.0  9207.0   9300.0\n  94.0  188.0  282.0  376.0  470.0  564.0     9118.0  9212.0  9306.0   9400.0\n  95.0  190.0  285.0  380.0  475.0  570.0     9215.0  9310.0  9405.0   9500.0\n  96.0  192.0  288.0  384.0  480.0  576.0  …  9312.0  9408.0  9504.0   9600.0\n  97.0  194.0  291.0  388.0  485.0  582.0     9409.0  9506.0  9603.0   9700.0\n  98.0  196.0  294.0  392.0  490.0  588.0     9506.0  9604.0  9702.0   9800.0\n  99.0  198.0  297.0  396.0  495.0  594.0     9603.0  9702.0  9801.0   9900.0\n 100.0  200.0  300.0  400.0  500.0  600.0     9700.0  9800.0  9900.0  10000.0"

    # test that no spurious visual lines are added when one element spans multiple lines
    v = fill!(Array{Any}(undef, 9), 0)
    v[1] = "look I'm wide! --- "^9
    r = replstr(v)
    @test startswith(r, "9-element Vector{Any}:\n  \"look I'm wide! ---")
    @test endswith(r, "look I'm wide! --- \"\n 0\n 0\n 0\n 0\n 0\n 0\n 0\n 0")

    # test vertical/diagonal ellipsis
    v = fill!(Array{Any}(undef, 50), 0)
    v[1] = "look I'm wide! --- "^9
    r = replstr(v)
    @test startswith(r, "50-element Vector{Any}:\n  \"look I'm wide! ---")
    @test endswith(r, "look I'm wide! --- \"\n 0\n 0\n 0\n 0\n 0\n 0\n 0\n 0\n 0\n ⋮\n 0\n 0\n 0\n 0\n 0\n 0\n 0\n 0\n 0")

    r = replstr([fill(0, 50) v])
    @test startswith(r, "50×2 Matrix{Any}:\n 0  …   \"look I'm wide! ---")
    @test endswith(r, "look I'm wide! --- \"\n 0     0\n 0     0\n 0     0\n 0     0\n 0  …  0\n 0     0\n 0     0\n 0     0\n 0     0\n ⋮  ⋱  \n 0     0\n 0     0\n 0     0\n 0     0\n 0  …  0\n 0     0\n 0     0\n 0     0\n 0     0")

    # issue #34659
    @test replstr(Int32[]) == "Int32[]"
    @test replstr([Int32[]]) == "1-element Vector{Vector{Int32}}:\n []"
    @test replstr(permutedims([Int32[], Int32[]])) == "1×2 Matrix{Vector{Int32}}:\n []  []"
    @test replstr(permutedims([Dict(), Dict()])) == "1×2 Matrix{Dict{Any, Any}}:\n Dict()  Dict()"
    @test replstr(permutedims([undef, undef])) == "1×2 Matrix{UndefInitializer}:\n UndefInitializer()  UndefInitializer()"
    @test replstr([zeros(3, 0), zeros(2, 0)]) == "2-element Vector{Matrix{Float64}}:\n 3×0 Matrix{Float64}\n 2×0 Matrix{Float64}"
end
@test replstr(Int32[]) == "Int32[]"
@test replstr([Int32[]]) == "1-element Vector{Vector{Int32}}:\n []"

# of a subarray

a = rand(5 * 10^2, 5 * 10^2) #5 x 5
s = view(a, 2:3, 2:3) # view
p = permutedims(s, [2, 1]) #permutation
@benchmark sum(p)

"""a = rand(5 * 10^2, 5 * 10^2)
27.711 ns (1 allocation: 16 bytes)
2.776234565017754

BenchmarkTools.Trial: 10000 samples with 996 evaluations.
 Range (min … max):  27.410 ns …  2.695 μs  ┊ GC (min … max): 0.00% … 97.07%
 Time  (median):     36.747 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   46.968 ns ± 45.473 ns  ┊ GC (mean ± σ):  1.31% ±  1.69%

  █▇▁▁▃▆▂▂▃▁▁▂▁▂▄▅▂▁▂▁▂▂▃▂▂▁▁▁▁▁    ▁▁▁                       ▂
  ██████████████████████████████████████▇▇▇▇▇▅▆▆▅▅▅▅▅▅▆▅▅▄▅▄▅ █
  27.4 ns      Histogram: log(frequency) by time       136 ns <

 Memory estimate: 16 bytes, allocs estimate: 1.

 5 * 10, 5 * 10
BenchmarkTools.Trial: 10000 samples with 995 evaluations.
 Range (min … max):  25.729 ns …  2.042 μs  ┊ GC (min … max): 0.00% … 97.26%
 Time  (median):     28.744 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   36.749 ns ± 38.906 ns  ┊ GC (mean ± σ):  1.47% ±  1.69%

   ▄█▅▁▃▅▄▂ ▁▂▂▁ ▁▁   ▁                                        ▁
  █████████████████████▇▇███▇███▇▇▇▆▇▆▆▆▆▅▅▅▅▅▅▅▃▄▃▄▅▅▅▅▆▅▅▅▅ █
  25.7 ns      Histogram: log(frequency) by time       115 ns <

 Memory estimate: 16 bytes, allocs estimate: 1.

"""
# ok
# of a non-strided subarray
a = reshape(1:60, 3, 4, 5)
s = view(a, :, [1, 2, 4], [1, 5])
c = convert(Array, s) #ok

for p in ([1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1])
    @test permutedims(s, p) == permutedims(c, p)
    @test PermutedDimsArray(s, p) == permutedims(c, p) #ok
end
@test_throws ArgumentError permutedims(a, (1, 1, 1))
@test_throws ArgumentError permutedims(s, (1, 1, 1))
@test_throws ArgumentError PermutedDimsArray(a, (1, 1, 1))
@test_throws ArgumentError PermutedDimsArray(s, (1, 1, 1))
cp = PermutedDimsArray(c, (3, 2, 1))

v = view(A, :, 3, 2:5) # bounds error
@test summary(v) == "2×4 view(::Array{Int16, 3}, :, 3, 2:5) with eltype Int16"
@test Base.showarg(io, v, false) === nothing
@test String(take!(io)) == "view(::Array{Int16, 3}, :, 3, 2:5)"

r = reshape(v, 4, 2) #reshape inconsistent
@test summary(r) == "4×2 reshape(view(::Array{Int16, 3}, :, 3, 2:5), 4, 2) with eltype Int16"
@test Base.showarg(io, r, false) === nothing
@test String(take!(io)) == "reshape(view(::Array{Int16, 3}, :, 3, 2:5), 4, 2)"

@testset "showarg" begin
    io = IOBuffer()

    A = reshape(Vector(Int16(1):Int16(2 * 3 * 5)), 2, 3, 5)
    @test summary(A) == "2×3×5 Array{Int16, 3}"

    v = view(A, :, 3, 2:5)
    @test summary(v) == "2×4 view(::Array{Int16, 3}, :, 3, 2:5) with eltype Int16"
    @test Base.showarg(io, v, false) === nothing
    @test String(take!(io)) == "view(::Array{Int16, 3}, :, 3, 2:5)"

    r = reshape(v, 4, 2)
    @test summary(r) == "4×2 reshape(view(::Array{Int16, 3}, :, 3, 2:5), 4, 2) with eltype Int16"
    @test Base.showarg(io, r, false) === nothing
    @test String(take!(io)) == "reshape(view(::Array{Int16, 3}, :, 3, 2:5), 4, 2)"

    p = PermutedDimsArray(r, (2, 1))
    @test summary(p) == "2×4 PermutedDimsArray(reshape(view(::Array{Int16, 3}, :, 3, 2:5), 4, 2), (2, 1)) with eltype Int16"
    @test Base.showarg(io, p, false) === nothing
    @test String(take!(io)) == "PermutedDimsArray(reshape(view(::Array{Int16, 3}, :, 3, 2:5), 4, 2), (2, 1))"

    p = reinterpret(reshape, Tuple{Float32,Float32}, [1.0f0 3.0f0; 2.0f0 4.0f0])
    @test summary(p) == "2-element reinterpret(reshape, Tuple{Float32, Float32}, ::Matrix{Float32}) with eltype Tuple{Float32, Float32}"
    @test Base.showarg(io, p, false) === nothing
    @test String(take!(io)) == "reinterpret(reshape, Tuple{Float32, Float32}, ::Matrix{Float32})"

    r = Base.IdentityUnitRange(2:2)
    B = @view ones(2)[r]
    Base.showarg(io, B, false)
    @test String(take!(io)) == "view(::Vector{Float64}, $(repr(r)))"
end

@testset "Methods" begin
    m = which(sin, (Float64,))
    io = IOBuffer()
    show(io, "text/html", m)
    s = String(take!(io))
    @test occursin(" in Base.Math ", s)
end

module AlsoExportsPair
Pair = 0
export Pair
end

module TestShowType
export TypeA
struct TypeA end
#using ..AlsoExportsPair
end

B[3, 1, 2] == A[1, 2, 3] # true

Base.@propagate_inbounds function (Base.getindex{typeof(T, N, AA, perm, iperm)}(A::PermutedDimsArray{T,N,AA,perm,iperm},  #T errors our
    I::Vararg{Int,N}))
    getindex(A.parent, ntuple(d --> I[iperm[d]], Val{N})...) #TODO:arrow is just -->  T is not defined
end
"""
Base.@propagate_inbounds _getindex(i, a, parent...) = (a[i], _getindex(i, b'...)...)
_getindex(i) = ()
"""

m = 1;
n = 10;
a = ones(n, m)
b = ones(m, n)

#--- permutation

parent = [a b']
perm = sqrt.(parent)

iperm = perm .^ 2

struct PermutedDimsArray{T,N,AA<:AbstractArray,perm,iperm} <: AbstractArray{T,N}
    parent::AA
    dims::NTuple{N,Int}

end
P = PermutedDimsArray{a,n,parent,perm,iperm}

@time res1 = dotop1(a, b)
"""  0.691503 seconds (1.10 M allocations: 56.898 MiB, 9.34% gc time, 99.99% compilation time)
10×10 BitMatrix
on rerun:
 0.000039 seconds (2 allocations: 144 bytes)
10×10 BitMatrix
"""

using BenchmarkTools
function dotop2(a, b)
    @benchmark a .== b
end
res2 = dotop2(a, b)

"""
BenchmarkTools.Trial: 10000 samples with 136 evaluations.
 Range (min … max):  708.824 ns … 98.607 μs  ┊ GC (min … max): 0.00% … 98.60%
 Time  (median):     963.235 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.143 μs ±  1.564 μs  ┊ GC (mean ± σ):  2.79% ±  2.19%
   █
  ▃█▅▄▆▆▃▂▃▂▂▂▂▃▃▄▄▄▄▃▃▃▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  709 ns          Histogram: frequency by time         2.78 μs <

 Memory estimate: 208 bytes, allocs estimate: 4.
"""
res3 = @btime a .== b # how to store it's values
"""
 rerun:

BenchmarkTools.Trial: 10000 samples with 138 evaluations.
 Range (min … max):  710.145 ns … 59.049 μs  ┊ GC (min … max): 0.00% … 96.91%
blue Time  (median):     841.304 ns              ┊ GC (median):    0.00%
green Time  (mean ± σ):     1.019 μs ±  1.329 μs  ┊ GC (mean ± σ):  2.98% ±  2.39%
  ▇█▆▅▅▅▅▃▄▄▄▃▃▃▃▃▃▃▃▃▂▂▁▁▁▁▁▁ ▁▁                              ▂
  ███████████████████████████████████▇█▇▇▇▇▆▇▆▇▆▆▅▅▆▅▄▆▄▃▅▄▄▄▄ █
  710 ns        Histogram: log(frequency) by time      2.47 μs <

 Memory estimate: 208 bytes, allocs estimate: 4.

where blue < green  (if there's a way to compare both distributions
including
the most frequency
blue line marker (median)
gren line marker (mean ± σ)  (alas, porvocu)
Notice: different runs yeilds hugely different discrepencies
(that may Hugely affect performance, on medium & long-run)
"""

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
module MVNTools

import POMDPs
import Base: rand, rand!
using AutoHashEquals

export MVN, apply_measurement!, apply_measurement, rand_elem

# Multivariate Normal Distribution
@auto_hash_equals type MVN <: POMDPs.AbstractDistribution
    mean::Vector{Float64}
    covariance::Array{Float64,2}
end
Base.length(d::MVN) = length(d.mean)

function apply_measurement!(d::MVN, index::Int, meas::Float64)
    @assert d.covariance[index,index] >= 0.0
    if d.covariance[index,index] == 0.0
        @assert meas == d.mean[index]
        return d
    end

    d.mean = d.mean + d.covariance[:,index]*(meas-d.mean[index])/d.covariance[index,index]

    d.covariance = d.covariance-d.covariance[:,index]*1/d.covariance[index,index]*d.covariance[index,:]

    d.covariance[:, index] = 0.0
    d.covariance[index, :] = 0.0

    return d
end

# the following version creates a copy
function apply_measurement(d::MVN, index::Int, meas::Float64)
    @assert d.covariance[index,index] >= 0.0
    @assert !isnan(meas)
    if d.covariance[index,index] == 0.0
        @assert meas == d.mean[index]
        return MVN(copy(d.mean), copy(d.covariance))
    end

    mean = d.mean + d.covariance[:,index]*(meas-d.mean[index])/d.covariance[index,index]

    covariance = d.covariance-d.covariance[:,index]*1/d.covariance[index,index]*d.covariance[index,:]

    covariance[:, index] = 0.0
    covariance[index, :] = 0.0

    # @assert all(mean.==mean) # this will check for NaNs

    return MVN(mean, covariance)
end

function rand!(rng::AbstractRNG, sample::Vector{Float64}, d::MVN; robust=false)
    # if robust, attempts to add a small number to the diagonal
    @assert length(sample) == length(d)

    nonzero = find(diag(d.covariance).!=0.0)
    nzcov = d.covariance[nonzero,nonzero]
    nzchol = zeros(size(nzcov))
    if robust
        try
            nzchol = chol(nzcov, Val{:L})
        catch e
            # Thu 19 Nov 2015 11:53:38 AM PST: couldn't remember why this is necessary
            nzchol = chol(nzcov+1e-5*minimum(diag(nzcov))*eye(length(nonzero)), Val{:L})
        end
    else
        nzchol = chol(nzcov, Val{:L})
    end
    delta = zeros(length(d))
    if rng == nothing
        delta[nonzero] = nzchol*randn(length(nonzero))
    else
        delta[nonzero] = nzchol*randn(rng, length(nonzero))
    end
    sample[:] = d.mean + delta
    return sample
end

function rand(rng::AbstractRNG, d::MVN)
    sample = Array(Float64, length(d))
    rand!(rng, sample, d)
end

rand(d::MVN) = rand(Base.GLOBAL_RNG, d)

# returns a single element of a random sample
function rand_elem(rng::AbstractRNG, d::MVN, elem::Int; robust=false)
    return d.mean[elem] + sqrt(d.covariance[elem, elem])*randn(rng)
end

#= this should be taken care of by auto hash equals
==(u::MVN, v::MVN) = u.mean==v.mean && u.covariance==v.covariance
hash(u::MVN, h::Int) = hash(u.mean, hash(u.covariance, h))
=#

end


text:
const ModelBoundaryConditions = NamedTuple{(:solution, :tendency, :pressure)}

function ModelBoundaryConditions(tracers, proposal_bcs::NamedTuple)
    solution_boundary_conditions = SolutionBoundaryConditions(tracers, proposal_bcs)

    model_boundary_conditions = (solution = solution_boundary_conditions,
                                 tendency = TendenciesBoundaryConditions(solution_boundary_conditions),
                                 pressure = PressureBoundaryConditions(solution_boundary_conditions.v))
    return model_boundary_conditions
end

ModelBoundaryConditions(tracers, model_boundary_conditions::ModelBoundaryConditions) =
    model_boundary_conditions


text:
mutable struct DESPOTConfig
  # Maximum depth of the search tree
  search_depth::Int64
  # Random-number seed
  main_seed::UInt32
  # Amount of CPU time used for search during each move. Does not include the
  # time taken to prune the tree and update the belief.
  time_per_move::Float64
  # Number of starting states (samples)
  n_particles::Int64
  # Regularization parameter
  pruning_constant::Float64
  # Parameter such that eta * width(root) is the target uncertainty at the
  # root of the search tree, used in determining when to terminate a trial.
  eta::Float64
  # Number of moves to simulate
  sim_len::Int64
  # Whether the initial upper bound is approximate or not. If approximate,
  # the solver allows initial lower bound > initial upper bound at a node.
  approximate_ubound::Bool
#   particle_weight_threshold::Float64
#   eff_particle_fraction::Float64
  tiny::Float64 # tiny number
  max_trials::Int64
  rand_max::Int64
  debug::UInt8

  # construct empty
  function DESPOTConfig()
   this = new()
   return this
  end
end


text:
using RvSpectraKitLearn

function make_B(p::Integer, q::Integer; frac_active::Float64 = 0.1, bias::Float64 = 0.0, sigma::Float64 = 1.0)
  B = zeros(p,q)
  for i in 1:p, j in 1:q
      if rand() < frac_active
        B[i,j] = bias+sigma*randn()
      end
  end
  return B
end

# Test sparse regression via  alg
n = 25
 p = 100
 q = 4
 frac_active = 0.1
 sigma_x = 0.1
 sigma_y = 0.1
 X = randn(n,p)
 B = make_B(p,q,frac_active=frac_active,bias=1.0,sigma=sigma_x)
 Y = X*B + sigma_y*randn(n,q)
 X_crossvalid = randn(n,p)
 Y_crossvalid = X_crossvalid*B + sigma_y*randn(n,q)
 println("B norms: ", length(find(B)), "  ", sumabs(B), " ", sumabs2(B))

(gamma_list, Z_list)  = fit_arp_slr_2d(X,Y,stepsize=1.04,max_itterations=1000, min_active=floor(Int64,p*q*frac_active*0.1))
(L0,L1,L2,rmsdelta) = compute_alr_path_stats(Z_list,target=B)
(chisq,chisq_refit,chisq_cv) = compute_alr_path_chisq!(Z_list,X,Y,X_crossvalid,Y_crossvalid)
print_alr_path(gamma_list, Z_list,X,Y)

#=
#Pkg.add("PyPlot")  # Run once if you don't already have PyPlot installed
using PyPlot

plt[:loglog](gamma_list,L0,"y-")
 plt[:loglog](gamma_list,L1/length(B),"r-")
 plt[:loglog](gamma_list,sqrt(L2),"g-")
 plt[:loglog](gamma_list,chisq_refit,"b-")
 plt[:loglog](gamma_list,chisq_refit,"b.")
 plt[:loglog](gamma_list,chisq_cv,"m-")
 plt[:loglog](gamma_list,chisq_cv,"m.")

plot(L0,chisq,"r-")
 plot(L0,chisq,"r.")
 plot(L0,chisq_refit,"b-")
 plot(L0,chisq_refit,"b.")
 plot(L0,chisq_cv,"m-")
 plot(L0,chisq_cv,"m.")
 plt[:loglog](L0,chisq_cv,"m.")

plot(L1,chisq,"r-")
 plot(L1,chisq,"r.")
 plot(L1,chisq_refit,"b-")
 plot(L1,chisq_refit,"b.")
 plot(L1,chisq_cv,"m-")
 plot(L1,chisq_cv,"m.")
 plt[:loglog](L1,chisq_cv,"m.")

plot(L0,L1./L0,"g-")
 plot(L0,L1./L0,"g.")

plot(L0,chisq,"r-")
 plot(L0,chisq,"r.")
 plot(L0,chisq_refit,"b-")
 plot(L0,chisq_refit,"b.")
 plot(L0,chisq_cv,"m-")
 plot(L0,chisq_cv,"m.")
 plt[:loglog](L0,chisq_cv,"m.")

=#

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
# ## Line plot 2d

#md # ![FILE_NAME.png](images/FILE_NAME.png)

x = 0:0.1:2pi
y1 = cos.(x)
y2 = sin.(x)

plot(x, y1, c="blue", linewidth=3)
plot!(x, y2, c="red", line=:dash)
title!("Trigonometric functions")
xlabel!("angle")
ylabel!("sin(x) and cos(x)")

## axis limit
plot!(xlims=(0,2pi), ylims=(-2, 2))

#-
## same plot
x = 0:0.1:2pi
y1 = cos.(x)
y2 = sin.(x)

plot(x, y1,
    c="blue",
    linewidth=3,
    title="Trigonometric functions",
    xlabel="angle",
    ylabel="sin(x) and cos(x)")
plot!(x, y2, c="red", line=:dash)

plot!(xlims=(0,2pi), ylims=(-2, 2))


text:
using MakieGallery, AbstractPlotting, GLMakie

push!(MakieGallery.plotting_backends, "GLMakie")

database = MakieGallery.load_database()
# THese examples download additional data - don't want to deal with that!
to_skip = ["WorldClim visualization", "Image on Geometry (Moon)", "Image on Geometry (Earth)"]
# we directly modify the database, which seems easiest for now
filter!(entry-> !(entry.title in to_skip), database)

ref_path = MakieGallery.download_reference(v"0.0.9")

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)
rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)

MakieGallery.record_examples(test_record_path)
MakieGallery.run_comparison(test_record_path, ref_path, tested_diff_path)


text:
using LifeTable
using Test
using DataFrames
using CSV
using GLM

# Builds a lifetable for Swedish males 2014 (based on data from Statistics Sweden).
mse14 = CSV.read("data/mse14.csv")
ltmse14 = periodlifetable(mse14, 1)

@test isapprox(ltmse14[:e][1], 80.3453; atol=1e-3)

# Performs a regression on Weibull parameters for survival for Swedish
# females dying of leukemia in age intervals from 30--34 to 80--84 years
# 1951--2012 (based on data from WHO Mortality Database).

leuk_msfs = Dict[]
fse5112tot3084 = CSV.read("data/fse5112tot3084.csv")
fse5112pop3084 = CSV.read("data/fse5112pop3084.csv")
fse5112leuk3084 = CSV.read("data/fse5112leuk3084.csv")

for i in 2:63
	dths = DataFrame(age = fse5112pop3084[1], p = fse5112pop3084[i], d = fse5112tot3084[i])
	lt = periodlifetable(dths, 2, false)
	leuk_prop = fse5112leuk3084[i] ./ fse5112tot3084[i]
	leuk_life = causelife(lt, leuk_prop)
	leuk_lt = periodlifetable(leuk_life, 2, false, "rate")
	leuk_msf = mortsurvfit(leuk_lt, fse5112tot3084[i], "weibull", "surv")
	push!(leuk_msfs, leuk_msf)
end

fit_leuk_msfs = mortsurvparamsfit(leuk_msfs)

@test isapprox(coef(fit_leuk_msfs["fit"])[2], -4.5; atol=1e-1)


text:
# The Numerical stability checks code in this file is based on the code from the SDDP.jl package,
# from the below mentioned commit and file.
# commit :8cd305188caffc50a1734913053fc81bba613778
# link to file :https://github.com/odow/SDDP.jl/blob/d353fe5a2903421e7fed6d609eb9377c35d715a1/src/print.jl#L190

mutable struct NumericalBounds
    min::Float64
    max::Float64
    min_index::Any
    max_index::Any
end

NumericalBounds() = NumericalBounds(Inf, -Inf, nothing, nothing)

set_min!(v::NumericalBounds, value::Real) = v.min = value
set_max!(v::NumericalBounds, value::Real) = v.max = value
set_min_index!(v::NumericalBounds, idx) = v.min_index = idx
set_max_index!(v::NumericalBounds, idx) = v.max_index = idx

mutable struct ConstraintBounds
    coefficient::NumericalBounds
    rhs::NumericalBounds
    function ConstraintBounds()
        return new(NumericalBounds(), NumericalBounds())
    end
end

function update_coefficient_bounds(
    v::ConstraintBounds,
    constraint::JuMP.ScalarConstraint,
    idx,
)
    update_numerical_bounds(v.coefficient, constraint.func, idx)
    return
end

function update_rhs_bounds(v::ConstraintBounds, constraint::JuMP.ScalarConstraint, idx)
    update_numerical_bounds(v.rhs, constraint.set, idx)
    return
end

mutable struct VariableBounds
    bounds::NumericalBounds
    function VariableBounds()
        return new(NumericalBounds())
    end
end

function update_variable_bounds(v::VariableBounds, variable::JuMP.VariableRef, idx)
    if JuMP.is_binary(variable)
        set_min!(v.bounds, 0.0)
        update_numerical_bounds(v.bounds, 1.0, idx)
    else
        if JuMP.has_lower_bound(variable)
            update_numerical_bounds(v.bounds, JuMP.lower_bound(variable), idx)
        end
        if JuMP.has_upper_bound(variable)
            update_numerical_bounds(v.bounds, JuMP.upper_bound(variable), idx)
        end
    end
    return
end

function update_numerical_bounds(v::NumericalBounds, value::Real, idx)
    if !isapprox(value, 0.0)
        if v.min > abs(value)
            set_min!(v, value)
            set_min_index!(v, idx)
        elseif v.max < abs(value)
            set_max!(v, value)
            set_max_index!(v, idx)
        end
    end
    return
end

function update_numerical_bounds(bonuds::NumericalBounds, func::JuMP.GenericAffExpr, idx)
    for coefficient in values(func.terms)
        update_numerical_bounds(bonuds, coefficient, idx)
    end
    return
end

function update_numerical_bounds(bonuds::NumericalBounds, func::MOI.LessThan, idx)
    return update_numerical_bounds(bonuds, func.upper, idx)
end

function update_numerical_bounds(bonuds::NumericalBounds, func::MOI.GreaterThan, idx)
    return update_numerical_bounds(bonuds, func.lower, idx)
end

function update_numerical_bounds(bonuds::NumericalBounds, func::MOI.EqualTo, idx)
    return update_numerical_bounds(bonuds, func.value, idx)
end

function update_numerical_bounds(bonuds::NumericalBounds, func::MOI.Interval, idx)
    update_numerical_bounds(bonuds, func.upper, idx)
    return update_numerical_bounds(bonuds, func.lower, idx)
end

# Default fallback for unsupported constraints.
update_numerical_bounds(::NumericalBounds, func, idx) = nothing

function get_constraint_numerical_bounds(model::OperationModel)
    if !is_built(model)
        error("Model not built, can't calculate constraint numerical bounds")
    end
    bounds = ConstraintBounds()
    for (const_key, constriant_array) in get_constraints(get_optimization_container(model))
        # TODO: handle this at compile and not at run time
        if isa(constriant_array, SparseAxisArray)
            for idx in eachindex(constriant_array)
                constriant_array[idx] == 0.0 && continue
                con_obj = JuMP.constraint_object(constriant_array[idx])
                update_coefficient_bounds(bounds, con_obj, (const_key, idx))
                update_rhs_bounds(bounds, con_obj, (const_key, idx))
            end
        else
            for idx in Iterators.product(constriant_array.axes...)
                !isassigned(constriant_array, idx...) && continue
                con_obj = JuMP.constraint_object(constriant_array[idx...])
                update_coefficient_bounds(bounds, con_obj, (const_key, idx))
                update_rhs_bounds(bounds, con_obj, (const_key, idx))
            end
        end
    end
    return bounds
end

function get_variable_numerical_bounds(model::OperationModel)
    if !is_built(model)
        error("Model not built, can't calculate variable numerical bounds")
    end
    bounds = VariableBounds()
    for (variable_key, variable_array) in get_variables(get_optimization_container(model))
        if isa(variable_array, SparseAxisArray)
            for idx in eachindex(variable_array)
                var = variable_array[idx]
                var == 0.0 && continue
                update_variable_bounds(bounds, var, (variable_key, idx))
            end
        else
            for idx in Iterators.product(variable_array.axes...)
                var = variable_array[idx...]
                update_variable_bounds(bounds, var, (variable_key, idx))
            end
        end
    end
    return bounds
end


text:
function helper_fn(z, attempt)
	single_data_name, test_fold = z
	proc_id = myid();

	helper_start_time = time()

	transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_transaction_counts.csv") ) )
	assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_assortments.csv") ) )
	folds  = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_folds.csv") ) )

	# M = size(assortments,1)

	# folds = collect(1:M)

	# train_transactions = transactions[ find(folds .!= 1), :]
	# test_transactions = transactions[ find(folds .== 1), :]

	train_assortments = assortments[ find(folds .!= test_fold), :]
	test_assortments = assortments[ find(folds .== test_fold), :]

	train_transaction_counts = transaction_counts[ find(folds .!= test_fold), :]
	test_transaction_counts = transaction_counts[ find(folds .== test_fold), :]

	num_train_transactions = sum(train_transaction_counts)
	num_test_transactions = sum(test_transaction_counts)

	num_train_assortments = size(train_assortments, 1)
	num_test_assortments = size(test_assortments, 1)

	N = 10

	time_limit = 1200

	perTol = 1e-6

	## RANKING-BASED MODEL ##

	ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking(N, train_assortments, train_transaction_counts, time_limit, perTol)
	ranking_K = length(ranking_lambda)

	test_ranking_predict, test_ranking_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, test_assortments)
	ranking_KL, ranking_noinf_count, ranking_KL_avg = tcm_evaluateKL(N, test_ranking_predict, test_transaction_counts)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"ranking", ",",
							ranking_loglik, ",",
							ranking_K, ",",
							perTol, ",",
							ranking_KL, ",",
							ranking_noinf_count, ",",
							ranking_KL_avg, ",",
							ranking_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(ranking_outcsvhandle, row_string);
		flush(ranking_outcsvhandle);
	end

	## LC-MNL MODEL ##

	K_list = [2, 3, 4, 5, 10, 15]
	numInitializations = 5

	for K in K_list
		LCMNL_u_val, LCMNL_class_probs, LCMNL_loglik, LCMNL_grand_elapsed_time = tcm_estimateMMNL(N, K, train_assortments, train_transaction_counts, numInitializations, time_limit, perTol)

		test_LCMNL_predict, test_LCMNL_predict_by_class = tcm_predictMMNL(N, K, LCMNL_class_probs, LCMNL_u_val, test_assortments)
		LCMNL_KL, LCMNL_noinf_count, LCMNL_KL_avg = tcm_evaluateKL(N, test_LCMNL_predict, test_transaction_counts)

		row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"LCMNL", ",",
							LCMNL_loglik, ",",
							K, ",",
							perTol, ",",
							numInitializations, ",",
							LCMNL_KL, ",",
							LCMNL_noinf_count, ",",
							LCMNL_KL_avg, ",",
							LCMNL_grand_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(MNL_outcsvhandle, row_string);
			flush(MNL_outcsvhandle);
		end
	end

	## LCMNL-CV MODEL ##

	numInitializations = 5
	perTol_CV = 5e-6;
	K_list = [2,3,4,5,10,15]
	nfolds = 4
	train_folds = folds[ find(folds .!= test_fold)]
	best_K, best_KL_cv, cv_elapsed_time = tcm_estimateMMNL_kcv_weighted_ass(N, K_list, nfolds, train_folds, train_assortments, train_transaction_counts, numInitializations, time_limit, perTol_CV)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"LCMNL-CV", ",",
							best_K, ",",
							perTol_CV, ",",
							numInitializations, ",",
							best_KL_cv, ",",
							cv_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(LCMNL_cv_outcsvhandle, row_string);
		flush(LCMNL_cv_outcsvhandle);
	end

	## MNL MODEL ##

	MNL_u_val, MNL_loglik, MNL_elapsed_time = tcm_estimateMNL(N, train_assortments, train_transaction_counts)
	test_MNL_predict = tcm_predictMNL(N, MNL_u_val, test_assortments)
	MNL_KL, MNL_noinf_count, MNL_KL_avg = tcm_evaluateKL(N, test_MNL_predict, test_transaction_counts)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"MNL", ",",
							MNL_loglik, ",",
							1, ",",
							0, ",",
							1, ",",
							MNL_KL, ",",
							MNL_noinf_count, ",",
							MNL_KL_avg, ",",
							MNL_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(MNL_outcsvhandle, row_string);
		flush(MNL_outcsvhandle);
	end

	## HALO-MNL MODEL ##

	HALOMNL_u_val, HALOMNL_alpha_val, HALOMNL_loglik, HALOMNL_elapsed_time = tcm_estimateHALOMNL(N, train_assortments, train_transaction_counts)

	test_HALOMNL_predict = tcm_predictHALOMNL(N, HALOMNL_u_val, HALOMNL_alpha_val, test_assortments)

	HALOMNL_KL, HALOMNL_noinf_count, HALOMNL_KL_avg = tcm_evaluateKL(N, test_HALOMNL_predict, test_transaction_counts)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"HALOMNL", ",",
							HALOMNL_loglik, ",",
							1, ",",
							0, ",",
							1, ",",
							HALOMNL_KL, ",",
							HALOMNL_noinf_count, ",",
							HALOMNL_KL_avg, ",",
							HALOMNL_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(MNL_outcsvhandle, row_string);
		flush(MNL_outcsvhandle);
	end

	## FOREST, COLD-START ##

	depth_list = [3,4,5,6,7]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_createBasicForest(N)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_loglik, forest_elapsed_time = tcm_estimateForest(N, train_assortments, train_transaction_counts, time_limit, perTol, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		test_forest_predict, test_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_KL, forest_noinf_count, forest_KL_avg = tcm_evaluateKL(N, test_forest_predict, test_transaction_counts)

		row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cold", ",",
							depth_limit, ",",
							forest_loglik, ",",
							K, ",",
							perTol, ",",
							forest_KL, ",",
							forest_noinf_count, ",",
							forest_KL_avg, ",",
							forest_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	## FOREST, WARM-START ##
	depth_list = [3,4,5,6,7]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_convertRankingToForest(N, orderings)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_loglik, forest_elapsed_time = tcm_estimateForest(N, train_assortments, train_transaction_counts, time_limit, perTol, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		test_forest_predict, test_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_KL, forest_noinf_count, forest_KL_avg = tcm_evaluateKL(N, test_forest_predict, test_transaction_counts)

		row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-ranking-warm", ",",
							depth_limit, ",",
							forest_loglik, ",",
							K, ",",
							perTol, ",",
							forest_KL, ",",
							forest_noinf_count, ",",
							forest_KL_avg, ",",
							forest_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	## FOREST-CV, WARM-START ##
	depth_list = [3,4,5,6,7]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_convertRankingToForest(N, orderings)

	nfolds = 4
	train_folds = folds[ find(folds .!= test_fold)]
	best_d, best_KL_avg_cv, elapsed_time_cv = tcm_estimateForest_kcv_weighted_ass(N, depth_list, nfolds, train_folds, train_assortments, train_transaction_counts, time_limit, perTol, initial_left, initial_right, initial_product, initial_isLeaf)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cv-ranking-warm", ",",
							best_d, ",",
							perTol, ",",
							best_KL_avg_cv, ",",
							elapsed_time_cv)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(forest_cv_outcsvhandle, row_string);
		flush(forest_cv_outcsvhandle);
	end

	## FOREST-CV, COLD-START ##
	depth_list = [3,4,5,6,7]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_createBasicForest(N)

	nfolds = 4
	train_folds = folds[ find(folds .!= test_fold)]
	best_d, best_KL_avg_cv, elapsed_time_cv = tcm_estimateForest_kcv_weighted_ass(N, depth_list, nfolds, train_folds, train_assortments, train_transaction_counts, time_limit, perTol, initial_left, initial_right, initial_product, initial_isLeaf)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cv-cold", ",",
							best_d, ",",
							perTol, ",",
							best_KL_avg_cv, ",",
							elapsed_time_cv)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(forest_cv_outcsvhandle, row_string);
		flush(forest_cv_outcsvhandle);
	end

	helper_end_time = time()

	return helper_end_time - helper_start_time;

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
prepare_homotopy(H::AbstractHomotopy, alg::APCA{true}) = homogenize(H)
prepare_homotopy(H::AbstractHomotopy, alg::APCA{false}) = H

"""
    prepare_startvalue(H, startvalue, ::AbstractPredictorCorrectorAlgorithm{true})

Embeds a vector into the projective space if necessary, i.e. if it's length is one less
than the number of variables of `H`. `H` is assumed to be homogenized. After the (eventual)
embedding the value is normalized.
"""
function prepare_startvalue(H::AbstractHomotopy, startvalue, ::APCA{true})
    N = nvariables(H)
    n = length(startvalue)
    if N == n
        normalize(startvalue)
    elseif N - 1 == n
        normalize!([1; startvalue])
    else
        error("A start value has length $n. Excepted length $N or $(N-1).")
    end
end

function prepare_startvalue(H::AbstractHomotopy, startvalue, ::APCA{false})
    N = nvariables(H)
    n = length(startvalue)
    if N == n
        return startvalue
    else
        return error("A start_value has length $n. Excepted length $N.")
    end
end


text:
requiredTypes = [ "Assignment" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "..", "types", reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
## TESTING VEC METHODS ###############################################

a = Vec(1.0,2.0,3.0)
b = Vec(4.0,6.0,8.0)

@testset "Geometry: Vecs" begin
	@test isapprox(a,a)
	@test false ≈ isapprox(a,b)
	@test isapprox(a+b ,Vec(5.0,8.0,11.0) )
	@test isapprox(b-a ,Vec(3.0,4.0,5.0) )
	@test isapprox(2*a, Vec(2.0,4.0,6.0))
	@test isapprox(a*b, 40)
	@test isapprox(cross(a,b) , Vec(-2,4,-2))
end
## TESTING POINTS METHODS ###############################################

a = Point(1.0, 2.0, 3.0)
b = Point(4.0, 6.0, 8.0)

@testset "Geometry: Point" begin
	@test isapprox(a,a) ≈ true
	@test isapprox(a,b) ≈ false
	@test isapprox((a * 2.), Point(2.0, 4.0, 6.0)) ≈ true
	@test isapprox((a + b), Point(5.0, 8.0, 11.0)) ≈ true
	#@test isapprox((b - a), Vec(3.0, 4.0, 5.0)) ≈ true #uncomment when adding Vec - Vec method
end

## TESTING TRANSFORMATION METHODS ###############################################

T1 = Transformation([1.0 2.0 3.0 4.0 ; 5.0 6.0 7.0 8.0 ; 9.0 9.0 8.0 7.0 ; 6.0 5.0 4.0 1.0],
                    [-3.75 2.75 -1 0 ; 4.375 -3.875 2.0 -0.5 ; 0.5 0.5 -1.0 1.0 ; -1.375 0.875 0.0 -0.5] )

T2 = Transformation( [3.0 5.0 2.0 4.0 ; 4.0 1.0 0.0 5.0 ; 6.0 3.0 2.0 0.0 ; 1.0 4.0 2.0 1.0],
                     [0.4 -0.2 0.2 -0.6 ; 2.9 -1.7 0.2 -3.1 ; -5.55 3.15 -0.4 6.45 ; -0.9 0.7 -0.2 1.1] )

T1_same = Transformation( [1.0 2.0 3.0 4.0 ; 5.0 6.0 7.0 8.0 ; 9.0 9.0 8.0 7.0 ; 6.0 5.0 4.0 1.0],
						  [-3.75 2.75 -1 0 ; 4.375 -3.875 2.0 -0.5 ; 0.5 0.5 -1.0 1.0 ; -1.375 0.875 0.0 -0.5] )

T_expected = Transformation( [33.0 32.0 16.0 18.0 ; 89.0 84.0 40.0 58.0 ; 118.0 106.0 48.0 88.0 ; 63.0 51.0 22.0 50.0],
                     		 [-1.45 1.45 -1.0 0.6 ; -13.95 11.95 -6.5 2.6 ; 25.525 -22.025 12.25 -5.2 ; 4.825 -4.325 2.5 -1.1] )

T1_homogeneous = Transformation( [1.0 2.0 3.0 4.0 ; 5.0 6.0 7.0 8.0 ; 9.0 9.0 8.0 7.0 ; 0.0 0.0 0.0 1.0],
                                 [-3.75 2.75 -1 0 ; 5.75 -4.75 2.0 1.0 ; -2.25 2.25 -1.0 -2.0 ; 0.0 0.0 0.0 1.0] )

V_expected = Vec(14.0, 38.0, 51.0)
P_expected = Point(18.0, 46.0, 58.0)
N_expected = Normal(-8.75, 7.75, -3.0)

@testset "Geometry: Transformation" begin

    @test is_consistent(T1)
    @test isapprox(T1, T1_same)

    T1_same.m[1,3] += 100.0
    T1_same.invm[2,1] += 100.0

    @test isapprox(T1, T1_same) ≈ false
    @test isapprox(T1, T1_same) ≈ false

    @test is_consistent(T2)
    @test is_consistent(T_expected)
    @test is_consistent(T1_homogeneous)

    @test isapprox(T_expected, T1*T2)
    @test isapprox(V_expected, T1_homogeneous*Vec(1.0, 2.0, 3.0))
    @test isapprox(P_expected, T1_homogeneous*Point(1.0, 2.0, 3.0))
    @test isapprox(N_expected, T1_homogeneous*Normal(3.0, 2.0, 4.0))

    T1_copy = inverse(T1)
    prod = T1 * T1_copy

    @test is_consistent(T1_copy)
    @test is_consistent(prod)
    @test isapprox(prod, Transformation() )

    translation_1 = translation(Vec(1.0, 2.0, 3.0))
    translation_2 = translation(Vec(4.0, 6.0, 8.0))
    prod = translation_1 * translation_2
    expected = translation(Vec(5.0, 8.0, 11.0))

    @test is_consistent(translation_1)
    @test is_consistent(translation_2)
    @test is_consistent(prod)
    @test isapprox(expected, prod)

    @test is_consistent(rotation_y(0.1f0))
    @test is_consistent(rotation_z(0.1f0))
    @test isapprox( (rotation_x( pi/2. ) * Vec(0.0, 1.0, 0.0)) , (Vec(0.0, 0.0, 1.0)) )
    @test isapprox( (rotation_y( pi/2. ) * Vec(0.0, 0.0, 1.0)) , (Vec(1.0, 0.0, 0.0)) )
    @test isapprox( (rotation_z( pi/2. ) * Vec(1.0, 0.0, 0.0)) , (Vec(0.0, 1.0, 0.0)) )

end

pcg = PCG()

@testset "Geometry: Orthonormal bases" begin

    for i in range(1, length=10^3)
        normal = Vec(pcg_randf(pcg), pcg_randf(pcg), pcg_randf(pcg))
        normal = normalize(normal)
        e1, e2, e3 = create_onb(normal)

        @test e3 ≈ normal

        @test isapprox(1.0f0 , squared_norm(e1), atol = 5f-6)
        @test isapprox(1.0f0 , squared_norm(e2), atol = 5f-6)
        @test isapprox(1.0f0 , squared_norm(e3), atol = 5f-6)

        @test isapprox(0.0f0 , e1 * e2, atol = 5f-6)
        @test isapprox(0.0f0 , e2 * e3, atol = 5f-6)
        @test isapprox(0.0f0 , e3 * e1, atol = 5f-6)
    end
end


text:
# ------------------------------------------------------------
# We can use a barrier function
using BenchmarkTools

function double_sum(data)
    total = 0
    for v in data
        total += 2 * v
    end
    return total
end

function double_sum_of_random_data(n)
    data = random_data(n)
    return double_sum(data)
end

@btime double_sum_of_random_data(100000);
@btime double_sum_of_random_data(100001);

#=
julia> @btime double_sum_of_random_data(100000);
  245.044 μs (2 allocations: 781.33 KiB)

julia> @btime double_sum_of_random_data(100001);
  180.454 μs (2 allocations: 781.39 KiB)

=#

# Another potential issue is the accumulator.

@code_warntype double_sum(rand(Int, 3));
@code_warntype double_sum(rand(Float64, 3));
# -- red text for the Float64 call against the `total` variable
#=
julia> @code_warntype double_sum(rand(Int, 3));
Variables
  #self#::Core.Compiler.Const(double_sum, false)
  data::Array{Int64,1}
  total::Int64
  @_4::Union{Nothing, Tuple{Int64,Int64}}
  v::Int64

Body::Int64
1 ─       (total = 0)
│   %2  = data::Array{Int64,1}
│         (@_4 = Base.iterate(%2))
│   %4  = (@_4 === nothing)::Bool
│   %5  = Base.not_int(%4)::Bool
└──       goto #4 if not %5
2 ┄ %7  = @_4::Tuple{Int64,Int64}::Tuple{Int64,Int64}
│         (v = Core.getfield(%7, 1))
│   %9  = Core.getfield(%7, 2)::Int64
│   %10 = total::Int64
│   %11 = (2 * v)::Int64
│         (total = %10 + %11)
│         (@_4 = Base.iterate(%2, %9))
│   %14 = (@_4 === nothing)::Bool
│   %15 = Base.not_int(%14)::Bool
└──       goto #4 if not %15
3 ─       goto #2
4 ┄       return total

julia> @code_warntype double_sum(rand(Float64, 3));
Variables
  #self#::Core.Compiler.Const(double_sum, false)
  data::Array{Float64,1}
  total::Union{Float64, Int64}
  @_4::Union{Nothing, Tuple{Float64,Int64}}
  v::Float64

Body::Union{Float64, Int64}
1 ─       (total = 0)
│   %2  = data::Array{Float64,1}
│         (@_4 = Base.iterate(%2))
│   %4  = (@_4 === nothing)::Bool
│   %5  = Base.not_int(%4)::Bool
└──       goto #4 if not %5
2 ┄ %7  = @_4::Tuple{Float64,Int64}::Tuple{Float64,Int64}
│         (v = Core.getfield(%7, 1))
│   %9  = Core.getfield(%7, 2)::Int64
│   %10 = total::Union{Float64, Int64}
│   %11 = (2 * v)::Float64
│         (total = %10 + %11)
│         (@_4 = Base.iterate(%2, %9))
│   %14 = (@_4 === nothing)::Bool
│   %15 = Base.not_int(%14)::Bool
└──       goto #4 if not %15
3 ─       goto #2
4 ┄       return total

=#

# There are better ways to write it.

# Fix 1 - use zero function to match type
function double_sum(data)
    total = zero(eltype(data))
    for v in data
        total += 2 * v
    end
    return total
end

# check again - no more red text in the output
@code_warntype double_sum(rand(Int, 3));
@code_warntype double_sum(rand(Float64, 3));

# Fix 2 - leverage type parameter, more flexible
function double_sum(data::AbstractVector{T}) where {T <: Number}
    total = zero(T)
    for v in data
        total += v
    end
    return total
end

# check again - no more red text in the output
@code_warntype double_sum(rand(Int, 3));
@code_warntype double_sum(rand(Float64, 3));

# -----
# using zero(s), one(s), and similar

#=
julia> zero(Int)
0

julia> zeros(Float64, 5)
5-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0

julia> zeros(Float64, 5)^C

julia> one(UInt8)
0x01

julia> ones(UInt8, 5)
5-element Array{UInt8,1}:
 0x01
 0x01
 0x01
 0x01
 0x01

julia> A = rand(3,4)
3×4 Array{Float64,2}:
 0.275294  0.333407   0.380679  0.680078
 0.795881  0.297336   0.478146  0.377958
 0.534937  0.0791554  0.916814  0.210469

julia> B = similar(A)
3×4 Array{Float64,2}:
 9.88131e-324  4.94066e-324  4.94066e-324  4.94066e-324
 9.88131e-324  4.94066e-324  4.94066e-324  0.0
 4.94066e-324  4.94066e-324  4.94066e-324  2.32833e-314

julia> zeros(axes(A))
3×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

=#


text:
#=
    Contents: Create, read, update and delete functions for server-side sessions that use Redis as a database.

    Notes:
    1. Key paths are transformed from (k1, k2, k3) to k1:k2:k3.
    2. Schema:
       - "sessions" => Set(session_id1, ...),      set of current valid session_ids.
       - "session:$(session_id):$(keypath)" => v,  key-value pairs for session_id.
       - "session:keypaths" => Set(keypaths...),   exploits the fact that session fields come from a common app-specific set.
    3. Example:
       - set!(con, session_id, k1, k2, k3, v) = set(con, "session:$session_id:$k1:$k2:k3", v)
    4. Maintaining key_paths allows the delete_session! function to work.
=#
using Redis

"""
Init session_id => session on the database and set the specified cookie to session_id.
Return: session_id
"""
function session_create!(con::RedisConnection, username::AbstractString, res::Response, cookiename::AbstractString)
    get_n_sessions(username) >= config[:session][:max_n_sessions] && return
    session_id = generate_session_id()
    sadd(con, "sessions", session_id)                     # Add session_id to the "sessions" Set
    set(con, "session:$session_id:username", username)    # Add username to the session data
    sadd(con, "session:keypaths", "username")             # Add "username"  to the "session:keypaths" Set
    if haskey(config[:session], :timeout)                 # Add "lastvisit" to the "session:keypaths" Set
	sadd(con, "session:keypaths", "lastreq")
	update_lastreq!(con, session_id)
    end
    add_sessionid_to_user!(username, session_id)
    write_to_cookie!(res, cookiename, session_id)
    session_id
end

"Delete session from database and set the specified cookie to an invalid state."
function session_delete!(con::RedisConnection, session_id::AbstractString)
    if session_id != ""
	username = get(con, "session:$session_id:username")
	remove_sessionid_from_user!(username, session_id)
	delete_session_from_database!(con, session_id)
    end
end

function delete_session_from_database!(con::RedisConnection, session_id::AbstractString)
    # Delete session fields
    kys = smembers(con, "session:keypaths")
    for k in kys
	del(con, "session:$session_id:$k")
    end

    # Delete the session from the "sessions" Set
    srem(con, "sessions", session_id)
end

function delete_all_sessions_from_database!(con::RedisConnection)
    session_ids = smembers(con, "sessions")
    for session_id in session_ids
	delete_session_from_database!(con, session_id)
    end
end

function session_is_valid(con::RedisConnection, session_id::AbstractString)
    result = sismember(con, "sessions", session_id)
    if haskey(config[:session], :timeout)
        dt_str = session_get(con, session_id, "lastreq")
	if dt_str == nothing
	    result = false
	else
	    dt = DateTime(dt_str)
	    if dt + Dates.Second(config[:session][:timeout]) < now()
		result = false    # Session has timed out
	    end
	end
    end
    result
end

function update_lastreq!(con::RedisConnection, session_id::AbstractString)
    set(con, "session:$session_id:lastreq", string(now()))
end

"""
Returns: key path defined by keys..., converted to Redis format.

Example 1:
    - INPUT:  keys_value... = k1, k2, k3, v
    - OUTPUT: key_path = "k1:k2:k3"

Example 2:
    - INPUT:  keys_value... = k1, k2, k3 (NO VALUE...has_value == false)
    - OUTPUT: key_path = "k1:k2:k3"
"""
function construct_keypath(has_value::Bool, keys_value...)
    nkeys = length(keys_value)
    if has_value
	nkeys -= 1
    end
    assert(nkeys >= 1)
    io    = IOBuffer()
    print(io, "$(keys_value[1])")
    if nkeys >= 2
	for i = 2:nkeys
	    print(io, ":$(keys_value[i])")
	end
    end
    takebuf_string(io)
end

function session_get(con::RedisConnection, session_id, keys...)
    k = construct_keypath(false, keys...)
    get(con, "session:$session_id:$k")
end

function session_set!(con::RedisConnection, session_id, keys_value...)
    k = construct_keypath(true, keys_value...)
    sadd(con, "session:keypaths", k)                                      # Add key path to the set of key paths
    set(con, "session:$session_id:$k", keys_value[length(keys_value)])    # Set keypath => value
end

# EOF


text:
export readgrid

if VERSION < v"0.4.0-dev+1419"
  stringtoint(s)   = int(s);
  stringtofloat(s) = float(s);
else
  stringtoint(s)   = parse(Int64,s);
  stringtofloat(s) = parse(Float64,s);
end

"""
**Summary:**
```
 xyz,vert,pattr,ele,eattr,nd=function readgrid(fname::String)
```

**Description:**

Reads a *1D,2D,3D* topology grid from the file `fname` and stores data in
intermediate structure consisting of

- `xyz  ::Array{Float64,2}` vertices positions
- `vert ::Array{Int64,2}`   vertex numbers
- `pattr::Array{Int64,2}`   vertex attributes
- `ele  ::Array{Int64,2}`   element connectivity
- `eattr::Array{Int64,2}`   element attributes
- `nd::Int64`               spatial dimension (<=topological dimension)

 xyz,vert,pattr,ele,eattr,nd

as read from node and ele files. Can handle meshes generated with *triangle*,
*tetgen*, and the corresponding *1D* equivalent. The type of mesh (intervals,
triangles, tetraeders etc.) should be known by the caller of readgrid, when
creating the proper derived objects, i.e. topologies.

"""

function readgrid(fname::String)

  # read vertex related stuff
  fid = open(fname * ".node","r")
  mesh_header = split(strip(readline(fid))) # read file header
  nvert       = stringtoint(mesh_header[1]) # number of vertices
  nd          = stringtoint(mesh_header[2]) # spatial dimension
  npattr      = stringtoint(mesh_header[3]) # number of point attributes
  nidv        = stringtoint(mesh_header[4]) # number of boundary markers

  # allocate
  xyz   = zeros(Float64,nvert,nd)
  vert  = zeros(Int64, nvert, 2 )
  pattr = zeros(Int64, nvert, nidv + npattr)

  # read vertex positions, enumeration, and attributes
  for i=1:nvert
    data         = split(strip(readline(fid)))
    for j=1:nd
      xyz[i,j]     = stringtofloat(data[1+j])
    end
    vert[i,1]    = i
    vert[i,2]    = stringtoint(data[1])
    for j=1:nidv+npattr
      pattr[i,j] = stringtoint(data[1+nd+j])
    end
  end
  close(fid)

  # read element related stuff
  fid = open(fname * ".ele","r")
  mesh_header = split(strip(readline(fid)))
  nele   = stringtoint(mesh_header[1]) # number of elements
  ntop   = stringtoint(mesh_header[2]) # number of entities per element
  nide   = stringtoint(mesh_header[3]) # number of element attributes

  # allocate
  ele    = zeros(Int64,nele,ntop+1)
  eattr  = zeros(Int64,nele,nide)

  # read element enumeration, connectivity and attributes
  for k=1:nele
    data = split(strip(readline(fid)))
    for i=1:(ntop+1)
      ele[k,i] = stringtoint(data[i])
    end
    for i=1:nide
      eattr[k,i] = stringtoint(data[ntop+1+i])
    end
  end
  close(fid)
  return xyz,vert,pattr,ele,eattr,nd

end

# see http://userpages.umbc.edu/~squire/UCD_Form.htm
# readable with ParaView
function readgrid_avs(fname::String)

  # read vertex related stuff
  fid = open(fname * ".inp","r")
  mesh_header = split(strip(readline(fid))) # read file header
  nvert       = stringtoint(mesh_header[1]) # number of vertices
  nele        = stringtoint(mesh_header[2]) # number of elements
  npattr      = stringtoint(mesh_header[3]) # number of point attributes
  nide        = stringtoint(mesh_header[4]) # number of element attributes
  nmodel      = stringtoint(mesh_header[5]) # model data (assume this is zero)

  nd   = 3 # spatial dimension
  nide = 1 # overwrite above
  ntop = 3

  # allocate
  xyz   = zeros(Float64,nvert,nd)
  vert  = zeros(Int64, nvert, 2 )

  # read vertex positions, enumeration, and attributes
  for i=1:nvert
    data         = split(strip(readline(fid)))
    for j=1:nd
      xyz[i,j]     = stringtofloat(data[1+j])
    end
    vert[i,1]    = i
    vert[i,2]    = stringtoint(data[1])
  end

  ele    = zeros(Int64,nele,4)
  eattr  = zeros(Int64,nele,nide) # ==1?

  for k=1:nele
    data = split(strip(readline(fid)))
    if data[3]=="tri"
      ele[k,1] = stringtoint(data[1])
      for i=1:3
        ele[k,1+i] = stringtoint(data[3+i])
      end
      for i=1:nide
        eattr[k,i] = stringtoint(data[1+i])
      end
    end
  end

  # read from INP file?
  pattr = zeros(Int64,nvert,1)

  close(fid)
  return xyz,vert,pattr,ele,eattr,nd

end

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
abstract type AbstractGermination end

@columns struct ThresholdGermination{Mo} <: AbstractGermination
    # Field              | Def  | Unit | Pror            | Limits      | Log  | Description
    germination_mass::Mo | 1e-5 | mol  | Gamma(2.0, 2.0) | [1e-10,5.0] | true | "Structural mass at germination"
end

"""
Check if germination has happened. Independent for each organ,
although this may not make sense. A curve could be better for this too.
"""
is_germinated(o, u) = is_germinated(germination_pars(o), o, u)
is_germinated(f::Nothing, o, u) = true
is_germinated(f::ThresholdGermination, o, u) = u.V > f.germination_mass


text:
cd(dirname(@__FILE__))

module Temp

using Eratosthenes # Bring in the top-level sim environment.
using Printf

# Create the hierarchy of objects used by the sim.
scenario = setup("scenarios/reaction-wheel.yaml")

# Run it with text progress).
simulate(scenario) do k, n
    if (k-1) % 100 == 0
        @printf("Progress: % 5.1f%% done.\n", round(100*(k-1)/(n-1), digits=1))
    end
    return true
end

end # Temp


text:
#!/usr/bin/env julia --project

using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using Distributions
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cv_d, cp_d, MSLP, grav, LH_v0
using CLIMAParameters.Atmos.SubgridScale: C_smag
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.BalanceLaws:
    vars_state,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ClimateMachine.BalanceLaws: boundary_state!
import ClimateMachine.Atmos: flux_second_order!

"""
    Initial Condition for Taylor-Green vortex (LES)

## References:
 - [Taylor1937](@cite)
 - [Rafei2018](@cite)
 - [Bull2014](@cite)
"""
function init_greenvortex!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    # Problem float-type
    FT = eltype(state)
    param_set = parameter_set(bl)

    # Unpack constant parameters
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v

    # Compute perturbed thermodynamic state:
    ρ = FT(1.178)      # density
    ρu = SVector(FT(0), FT(0), FT(0))  # momentum
    #State (prognostic) variable assignment
    e_pot = FT(0)# potential energy
    Pinf = 101325
    Uzero = FT(100)
    p = Pinf + (ρ * Uzero^2 / 16) * (2 + cos(2 * z)) * (cos(2 * x) + cos(2 * y))
    u = Uzero * sin(x) * cos(y) * cos(z)
    v = -Uzero * cos(x) * sin(y) * cos(z)
    e_kin = 0.5 * (u^2 + v^2)
    T = p / (ρ * R_gas)
    e_int = internal_energy(param_set, T, PhasePartition(FT(0)))
    ρe_tot = ρ * (e_kin + e_pot + e_int)
    # Assign State Variables
    state.ρ = ρ
    state.ρu = SVector(FT(ρ * u), FT(ρ * v), FT(0))
    state.energy.ρe = ρe_tot
end

# Set up AtmosLESConfiguration for the experiment
function config_greenvortex(
    ::Type{FT},
    N,
    (xmin, xmax, ymin, ymax, zmin, zmax),
    resolution,
) where {FT}

    _C_smag = FT(C_smag(param_set))
    physics = AtmosPhysics{FT}(
        param_set;                          # Parameter set corresponding to earth parameters
        ref_state = NoReferenceState(),
        turbulence = Vreman(_C_smag),       # Turbulence closure model
        moisture = DryModel(),
    )

    model = AtmosModel{FT}(
        AtmosLESConfigType,                 # Flow in a box, requires the AtmosLESConfigType
        physics;                            # Atmos physics
        init_state_prognostic = init_greenvortex!,
        orientation = NoOrientation(),
        source = (),
    )

    config = ClimateMachine.AtmosLESConfiguration(
        "GreenVortex",          # Problem title [String]
        N,                      # Polynomial order [Int]
        resolution,             # (Δx, Δy, Δz) effective resolution [m]
        xmax,                   # Domain maximum size [m]
        ymax,                   # Domain maximum size [m]
        zmax,                   # Domain maximum size [m]
        param_set,              # Parameter set.
        init_greenvortex!,      # Function specifying initial condition
        boundary = ((0, 0), (0, 0), (0, 0)),
        periodicity = (true, true, true),
        xmin = xmin,
        ymin = ymin,
        zmin = zmin,
        model = model,                  # Model type
    )
    return config
end

# Define the diagnostics configuration for this experiment
function config_diagnostics(
    driver_config,
    (xmin, xmax, ymin, ymax, zmin, zmax),
    resolution,
    tnor,
    titer,
    snor,
)
    ts_dgngrp = setup_atmos_turbulence_stats(
        AtmosLESConfigType(),
        "360steps",
        driver_config.name,
        tnor,
        titer,
    )

    boundaries = [
        xmin ymin zmin
        xmax ymax zmax
    ]
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )
    ds_dgngrp = setup_atmos_spectra_diagnostics(
        AtmosLESConfigType(),
        "0.06ssecs",
        driver_config.name,
        interpol = interpol,
        snor,
    )
    me_dgngrp = setup_atmos_mass_energy_loss(
        AtmosLESConfigType(),
        "0.02ssecs",
        driver_config.name,
    )
    return ClimateMachine.DiagnosticsConfiguration([
        ts_dgngrp,
        ds_dgngrp,
        me_dgngrp,
    ],)
end

# Entry point
function main()
    FT = Float64
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Ncellsx = 64
    Ncellsy = 64
    Ncellsz = 64
    Δx = FT(2 * pi / Ncellsx)
    Δy = Δx
    Δz = Δx
    resolution = (Δx, Δy, Δz)
    xmin = FT(-pi)
    xmax = FT(pi)
    ymin = FT(-pi)
    ymax = FT(pi)
    zmin = FT(-pi)
    zmax = FT(pi)
    # Simulation time
    t0 = FT(0)
    timeend = FT(0.1)
    CFL = FT(1.8)

    driver_config = config_greenvortex(
        FT,
        N,
        (xmin, xmax, ymin, ymax, zmin, zmax),
        resolution,
    )

    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        init_on_cpu = true,
        Courant_number = CFL,
    )

    tnor = FT(100)
    titer = FT(0.01)
    snor = FT(10000.0)
    dgn_config = config_diagnostics(
        driver_config,
        (xmin, xmax, ymin, ymax, zmin, zmax),
        resolution,
        tnor,
        titer,
        snor,
    )

    check_cons = (
        ClimateMachine.ConservationCheck("ρ", "100steps", FT(0.0001)),
        ClimateMachine.ConservationCheck("energy.ρe", "100steps", FT(0.0025)),
    )

    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        check_cons = check_cons,
        check_euclidean_distance = true,
    )

end

main()


text:
module QJuliaComms

#list of error codes (opempi3.1.2)

const QJULIA_MPI_SUCCESS                   =0
const QJULIA_MPI_ERR_BUFFER                =1
const QJULIA_MPI_ERR_COUNT                 =2
const QJULIA_MPI_ERR_TYPE                  =3
const QJULIA_MPI_ERR_TAG                   =4
const QJULIA_MPI_ERR_COMM                  =5
const QJULIA_MPI_ERR_RANK                  =6
const QJULIA_MPI_ERR_REQUEST               =7
const QJULIA_MPI_ERR_ROOT                  =8
const QJULIA_MPI_ERR_GROUP                 =9
const QJULIA_MPI_ERR_OP                    =10
const QJULIA_MPI_ERR_TOPOLOGY              =11
const QJULIA_MPI_ERR_DIMS                  =12
const QJULIA_MPI_ERR_ARG                   =13
const QJULIA_MPI_ERR_UNKNOWN               =14
const QJULIA_MPI_ERR_TRUNCATE              =15
const QJULIA_MPI_ERR_OTHER                 =16
const QJULIA_MPI_ERR_INTERN                =17
const QJULIA_MPI_ERR_IN_STATUS             =18
const QJULIA_MPI_ERR_PENDING               =19
const QJULIA_MPI_ERR_ACCESS                =20
const QJULIA_MPI_ERR_AMODE                 =21
const QJULIA_MPI_ERR_ASSERT                =22
const QJULIA_MPI_ERR_BAD_FILE              =23
const QJULIA_MPI_ERR_BASE                  =24
const QJULIA_MPI_ERR_CONVERSION            =25
const QJULIA_MPI_ERR_DISP                  =26
const QJULIA_MPI_ERR_DUP_DATAREP           =27
const QJULIA_MPI_ERR_FILE_EXISTS           =28
const QJULIA_MPI_ERR_FILE_IN_USE           =29
const QJULIA_MPI_ERR_FILE                  =30
const QJULIA_MPI_ERR_INFO_KEY              =31
const QJULIA_MPI_ERR_INFO_NOKEY            =32
const QJULIA_MPI_ERR_INFO_VALUE            =33
const QJULIA_MPI_ERR_INFO                  =34
const QJULIA_MPI_ERR_IO                    =35
const QJULIA_MPI_ERR_KEYVAL                =36
const QJULIA_MPI_ERR_LOCKTYPE              =37
const QJULIA_MPI_ERR_NAME                  =38
const QJULIA_MPI_ERR_NO_MEM                =39
const QJULIA_MPI_ERR_NOT_SAME              =40
const QJULIA_MPI_ERR_NO_SPACE              =41
const QJULIA_MPI_ERR_NO_SUCH_FILE          =42
const QJULIA_MPI_ERR_PORT                  =43
const QJULIA_MPI_ERR_QUOTA                 =44
const QJULIA_MPI_ERR_READ_ONLY             =45
const QJULIA_MPI_ERR_RMA_CONFLICT          =46
const QJULIA_MPI_ERR_RMA_SYNC              =47
const QJULIA_MPI_ERR_SERVICE               =48
const QJULIA_MPI_ERR_SIZE                  =49
const QJULIA_MPI_ERR_SPAWN                 =50
const QJULIA_MPI_ERR_UNSUPPORTED_DATAREP   =51
const QJULIA_MPI_ERR_UNSUPPORTED_OPERATION =52
const QJULIA_MPI_ERR_WIN                   =53
const QJULIA_MPI_T_ERR_MEMORY              =54
const QJULIA_MPI_T_ERR_NOT_INITIALIZED     =55
const QJULIA_MPI_T_ERR_CANNOT_INIT         =56
const QJULIA_MPI_T_ERR_INVALID_INDEX       =57
const QJULIA_MPI_T_ERR_INVALID_ITEM        =58
const QJULIA_MPI_T_ERR_INVALID_HANDLE      =59
const QJULIA_MPI_T_ERR_OUT_OF_HANDLES      =60
const QJULIA_MPI_T_ERR_OUT_OF_SESSIONS     =61
const QJULIA_MPI_T_ERR_INVALID_SESSION     =62
const QJULIA_MPI_T_ERR_CVAR_SET_NOT_NOW    =63
const QJULIA_MPI_T_ERR_CVAR_SET_NEVER      =64
const QJULIA_MPI_T_ERR_PVAR_NO_STARTSTOP   =65
const QJULIA_MPI_T_ERR_PVAR_NO_WRITE       =66
const QJULIA_MPI_T_ERR_PVAR_NO_ATOMIC      =67
const QJULIA_MPI_ERR_RMA_RANGE             =68
const QJULIA_MPI_ERR_RMA_ATTACH            =69
const QJULIA_MPI_ERR_RMA_FLAVOR            =70
const QJULIA_MPI_ERR_RMA_SHARED            =71
const QJULIA_MPI_T_ERR_INVALID             =72
const QJULIA_MPI_T_ERR_INVALID_NAME        =73

#MPI communicator handles:
const QJULIA_MPI_COMM_WORLD = 0
const QJULIA_MPI_COMM_SELF  = 1
const QJULIA_MPI_COMM_NULL  = 2

mutable struct Comm
    val::Cint
    Comm(val::Integer) = new(val)
end

const COMM_NULL  = Comm(QJULIA_MPI_COMM_NULL)
const COMM_SELF  = Comm(QJULIA_MPI_COMM_SELF)
const COMM_WORLD = Comm(QJULIA_MPI_COMM_WORLD)

##### MPI & other communication calls ######
   MPI_init_qj(argc, argv)  = ccall((:MPI_Init, "libmpi"), Cvoid, (Int32, Ptr{Ptr{Cchar}}), argc, argv)

   MPI_initialized_qj(flag) = ccall((:MPI_Initialized, "libmpi"), Cint, (Ptr{Cint}, ), flag)

   MPI_rank_qj(comm, rank)  = ccall((:MPI_Comm_rank, "libmpi"), Cint, (Ref{Cint}, Ptr{Cint}), comm, rank)

   MPI_size_qj(comm, size)  = ccall((:MPI_Comm_size, "libmpi"), Cint, (Ref{Cint}, Ptr{Cint}), comm, size)

   MPI_finalize_qj()        = ccall((:MPI_Finalize, "libmpi"), Cvoid, (), )

##### High-level wrappers

function MPI_Initialized_qj()
   flag  = Ref{Cint}()
   error_code = MPI_initialized_qj(flag)
   if error_code != QJULIA_MPI_SUCCESS
     println("WARNING: MPI_Initialized_qj returned MPI error code ", error_code)
   end
   flag[] != 0
end

function MPI_Comm_rank_qj(comm::Comm)
   rank = Ref{Cint}()
   error_code = MPI_rank_qj(comm.val, rank)
   if error_code != QJULIA_MPI_SUCCESS
     println("WARNING: MPI_Comm_rank_qj returned MPI error code ", error_code)
   end
   Int(rank[])
end

function MPI_Comm_size_qj(comm::Comm)
   size = Ref{Cint}()
   error_code = MPI_size_qj(comm.val, size)
   if error_code != QJULIA_MPI_SUCCESS
     println("WARNING: MPI_Comm_size_qj returned MPI error code ", error_code)
   end
   Int(size[])
end

end #QJuliaComms

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
@testset "rank_correlation" begin
    for i in 1:1000
        len = rand(10:80)
        x = rand(1:100.0,len)
        y = rand(1:100.0,len)
        @test length(NighttimeLights.rank_correlation_test(x, y)) == 1
    end
end