text:
using SimpleTranslations
using Test

@testset "SimpleTranslations.jl" begin
    @test_throws ErrorException loadmsgs("notfound.ini")
    conf = loadmsgs("test.ini")
    @test get_msg(conf, "hi") == "Hi"
    set_language!(conf, "es")
    @test get_msg(conf, "hi") == "Hola a todos"
    loadmsgs("test_nostrict1.ini")
    @test_throws ErrorException loadmsgs("test_nostrict1.ini", strict_mode=true)
end

@testset "Usage with files" begin
    fname = joinpath(dirname(pathof(SimpleTranslations)), "..", "test", "test.ini")
    msgs = loadmsgs(fname)
    @test get_msg(msgs, "hi")=="Hi"
end

@testset "Checking strict" begin
    test_not_default = """
default = en

[en]
person=Person

[eo]
person=Homo
"""
    msgs = loadmsgs(IOBuffer(test_not_default), strict_mode=true)
    @test Set(["en", "eo"]) == Set(keys(msgs._msgs))
    @test Set(["en", "eo"]) == Set(msgs._languages)
end

@testset "Checking strict missing language" begin
    test_not_default = """
default = en
language=en,eo,es

[en]
person=Person

[eo]
person=Homo
"""
    @test_throws FileMessagesException msgs = loadmsgs(IOBuffer(test_not_default), strict_mode=true)
end

@testset "Checking strict missing language" begin
    test_not_default = """
default = en
language=en,eo,fr

[en]
person=Person
people=People

[eo]
person=Homo
people=Homoj

[fr]
person=Personne
"""
    @test_throws FileMessagesException msgs = loadmsgs(IOBuffer(test_not_default), strict_mode=true)
end

@testset "Global State" begin
    test_not_default = """
default = en
languages=en,eo,fr

[en]
person=Person
people=People

[eo]
person=Homo
people=Homoj

[fr]
person=Personne
people=Gens
"""
    loadmsgs!(IOBuffer(test_not_default), strict_mode=true)
    @test get_msg("person") == "Person"
    set_language!("fr")
    @test get_msg("person") == "Personne"
    set_language!("eo")
    @test get_msg("person") == "Homo"
end

@testset "Setting language" begin
    set_language!("en")
    fname = "test_strict.ini"
    loadmsgs!(fname, strict_mode=true)
    @test get_msg("hi") == "Hi"
    set_language!("es")
    loadmsgs!(fname, strict_mode=true)
    @test get_msg("hi") == "Hola a todos"
end

@testset "Setting initial language unknown" begin
    reset_msgs!()
    set_language!("eo")
    fname = "test_strict.ini"
    loadmsgs!(fname, strict_mode=false)
    @test get_msg("hi") == "Hi"
    reset_msgs!()
    set_language!("eo")
    @test_throws ErrorException loadmsgs!(fname, strict_mode=true)
end


text:
using JLD2
using Infiltrator
using LightGraphs
using Random
using JSON
using BenchmarkTools
using MultiAgentPathFinding

using AerialGroundPathFinding

const BASE_SEED = 6789
const MANHATTAN_NODEFILE = "../data/manhattan_node_attribs.json"
const MANHATTAN_WEIGHTFILE = "../data/manhattan_sparse_wts.jld2"

# Only caring about distance for now
const ALPHA_WEIGHT_DISTANCE = 1.0
const ECBS_WEIGHT = 1.5
# const THRESHOLD_CONFLICTS = 100

const NCARS = parse(Int64, ARGS[1])
const NDRONES = parse(Int64, ARGS[2])
const OUTFILEDIR = ARGS[3]
const CAPACITY = parse(Int64, ARGS[4])
const TRIALS = parse(Int64, ARGS[5])
const THRESHOLD_CONFLICTS = parse(Int64, ARGS[6])

function main()

    @load MANHATTAN_WEIGHTFILE manhattan_sparse_wts

    graph = LightGraphs.SimpleDiGraph(manhattan_sparse_wts)
    location_list = get_location2D_list(MANHATTAN_NODEFILE)

    for t = 1:TRIALS+1

        rng = MersenneTwister(BASE_SEED+t)

        ground_task_list = get_tasks_with_valid_path(graph, manhattan_sparse_wts, location_list, NCARS, rng)
        aerial_task_list = get_tasks_with_valid_path(graph, manhattan_sparse_wts, location_list, NDRONES, rng)

        env = CoordinatedMAPFEnv(ground_task_list=ground_task_list, aerial_task_list=aerial_task_list, road_graph=graph, road_graph_wts=manhattan_sparse_wts, location_list=location_list, alpha_weight_distance=ALPHA_WEIGHT_DISTANCE, car_capacity=CAPACITY, threshold_global_conflicts=THRESHOLD_CONFLICTS)

        time_comps = Float64[]

        t1 = @elapsed dir_drone_paths = compute_independent_paths(env, env.aerial_task_list)
        t2 = @elapsed dir_car_paths = compute_independent_paths(env, env.ground_task_list)
        push!(time_comps, t1+t2)

        t3 = @elapsed augment_road_graph_with_aerial_paths!(env, dir_drone_paths)
        push!(time_comps, t3)

        initial_gmapf_states = [GroundMAPFState(task.origin, false) for task in env.ground_task_list]

        t4 = @elapsed gmapf_solver = ECBSSolver{GroundMAPFState,GroundMAPFAction,Float64,SumOfCosts,GroundMAPFConflict,GroundMAPFConstraint,CoordinatedMAPFEnv}(env=env, weight=ECBS_WEIGHT)
        push!(time_comps, t4)

        ground_cbs_soln = PlanResult[]
        try
            t5 = @elapsed ground_cbs_soln = search!(gmapf_solver, initial_gmapf_states)
            push!(time_comps, t5)
        catch y
            if isa(y, DomainError)
                println("Too many conflicts in stage 1 of trial $(t-1)!")
                continue
            elseif  isa(y, BoundsError)
                println("ECBS ran out of nodes in stage 1 of trial $(t-1)")
                continue
            else
                throw(y)
            end
        end

        ground_conflicts = env.num_global_conflicts

        env.num_global_conflicts = 0

        t6 = @elapsed update_ground_paths_with_ground_mapf_result!(env, ground_cbs_soln)
        push!(time_comps, t6)
        stage1_time = sum(time_comps)

        t7 = @elapsed set_ground_transit_graph!(env, [ddp.total_dist for ddp in dir_drone_paths])
        push!(time_comps, t7)

        initial_amapf_states = [AerialMAPFState(idx=1, ground_transit_idx=sg[1]) for (i, sg) in enumerate(env.gtg_drone_start_goal_idxs)]

        t8 = @elapsed aerial_mapf_solver = ECBSSolver{AerialMAPFState,GroundTransitAction,Float64,SumOfCosts,GroundTransitConflict,GroundTransitConstraint,CoordinatedMAPFEnv}(env=env, weight=ECBS_WEIGHT)
        push!(time_comps, t8)

        aerial_cbs_soln = PlanResult[]
        try
            t9 = @elapsed aerial_cbs_soln = search!(aerial_mapf_solver, initial_amapf_states)
            push!(time_comps, t9)
        catch y
            if isa(y, DomainError)
                println("Too many conflicts in stage 2 of trial $(t-1)!")
                continue
            elseif  isa(y, BoundsError)
                println("ECBS ran out of nodes in stage 2 of trial $(t-1)")
                continue
            else
                throw(y)
            end
        end

        aerial_conflicts = env.num_global_conflicts

        t10 = @elapsed update_aerial_ground_paths_cbs_solution!(env, aerial_cbs_soln, dir_drone_paths)
        push!(time_comps, t10)

        not_coord_ids = count_cars_not_coordinating(env.ground_paths)
        for ncid in not_coord_ids
            env.ground_paths[ncid] = dir_car_paths[ncid]
        end

        total_dist = compute_total_cost(env, env.aerial_paths, env.ground_paths)

        total_time = sum(time_comps)
        stage2_time = total_time - stage1_time

        res_dict = Dict("compute_time" => total_time, "total_path_cost" => total_dist, "ground_conflicts" => ground_conflicts, "aerial_conflicts" => aerial_conflicts, "stage1_time" => stage1_time, "stage2_time" => stage2_time)

        if t > 1
            fname = string(OUTFILEDIR, "_Cars_", NCARS, "_Drones_", NDRONES, "_Capacity_", CAPACITY, "_Trial_", (t-1), ".json")
            println(fname)
            open(fname, "w") do f
                JSON.print(f, res_dict, 2)
            end
        end
    end
end

main()

ext: jl


text:
##### Beginning of file

Test.@test( Base.VERSION >= VersionNumber("1.0") )

Test.@test( RepoSnapshots.version() > VersionNumber(0) )

Test.@test(
    RepoSnapshots.version() ==
        RepoSnapshots.version(RepoSnapshots)
    )

Test.@test(
    RepoSnapshots.version() ==
        RepoSnapshots.version(first(methods(RepoSnapshots.eval)))
    )

Test.@test(
    RepoSnapshots.version() ==
        RepoSnapshots.version(RepoSnapshots.eval)
    )

Test.@test(
    RepoSnapshots.version() ==
        RepoSnapshots.version(RepoSnapshots.eval, (Any,))
    )

Test.@test( RepoSnapshots.version(TestModuleA) == VersionNumber("1.2.3") )

Test.@test( RepoSnapshots.version(TestModuleB) == VersionNumber("4.5.6") )

Test.@test_throws(
    ErrorException,
    RepoSnapshots._TomlFile(joinpath(mktempdir(),"1","2","3","4")),
    )

##### End of file


text:
function transformtest1()
  # a line of segments and a plane with nodes, holes
  # using all the units
  n1 = Node(1,0,0)
  n2 = Node(2,0,0)
  n3 = Node(0.003,0,0)
  n4 = Node(400,0,0)
  n5 = Node(5000,0,0)
  n6  = Node(6e6,0,0)
  n7 = Node(7/2.54e-2,0,0)
  n8 = Node(8/2.54e-5,0,0)
  g1 = Group([
    n1,
    Units("m"),
    n2,
    Units("km"),
    n3,
    Units("cm"),
    n4,
    Default(w=75, h=25),
    Units("mm"),
    n5,
    UniformPlane(x1=1000, y1=1000, z1=-1000,
                 x2=1000, y2=9000, z2=-1000,
                 x3=9000, y3=9000, z3=-1000,
                 thick = 200,
                 seg1 = 20,
                 seg2 = 20,
                 nodes = [Node(2000,2000,-1000),
                          Node(3000,3000,-1000)],
                 holes = [Point(4000,4000,-1000),
                          Rect(5000,5000,-1000,6000,6000,-1000),
                          Circle(7000,7000,-1000,500)]),
    Units("um"),
    n6,
    Units("in"),
    n7,
    Units("mils"),
    n8,
    Segment(n1,n2),
    Segment(n2,n3),
    Segment(n3,n4),
    Segment(n4,n5),
    Segment(n5,n6),
    Segment(n6,n7),
    Segment(n7,n8)
    ])
  verified =
  """
  N_1 x=1.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  .units m
  N_2 x=2.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  .units km
  N_3 x=3.000000000e-03 y=0.000000000e+00 z=0.000000000e+00
  .units cm
  N_4 x=4.000000000e+02 y=0.000000000e+00 z=0.000000000e+00
  .default
  +  w=7.500000000e+01 h=2.500000000e+01
  .units mm
  N_5 x=5.000000000e+03 y=0.000000000e+00 z=0.000000000e+00
  G_6
  + x1=1.000000000e+03 y1=1.000000000e+03 z1=-1.000000000e+03
  + x2=1.000000000e+03 y2=9.000000000e+03 z2=-1.000000000e+03
  + x3=9.000000000e+03 y3=9.000000000e+03 z3=-1.000000000e+03
  + thick=2.000000000e+02 seg1=20 seg2=20
  + N_7 (2.000000000e+03,2.000000000e+03,-1.000000000e+03)
  + N_8 (3.000000000e+03,3.000000000e+03,-1.000000000e+03)
  + hole point (4.000000000e+03, 4.000000000e+03, -1.000000000e+03)
  + hole rect (5.000000000e+03, 5.000000000e+03, -1.000000000e+03, 6.000000000e+03, 6.000000000e+03, -1.000000000e+03)
  + hole circle (7.000000000e+03, 7.000000000e+03, -1.000000000e+03, 5.000000000e+02)
  .units um
  N_9 x=6.000000000e+06 y=0.000000000e+00 z=0.000000000e+00
  .units in
  N_10 x=2.755905512e+02 y=0.000000000e+00 z=0.000000000e+00
  .units mils
  N_11 x=3.149606299e+05 y=0.000000000e+00 z=0.000000000e+00
  E_12 N_1 N_2
  E_13 N_2 N_3
  E_14 N_3 N_4
  E_15 N_4 N_5
  E_16 N_5 N_9
  E_17 N_9 N_10
  E_18 N_10 N_11
  """
  testelement(g1,verified)
  g2 = transform(g1,txyz(10,0,0))
  verified =
  """
  N_1 x=1.100000000e+01 y=0.000000000e+00 z=0.000000000e+00
  .units m
  N_2 x=1.200000000e+01 y=0.000000000e+00 z=0.000000000e+00
  .units km
  N_3 x=1.300000000e-02 y=0.000000000e+00 z=0.000000000e+00
  .units cm
  N_4 x=1.400000000e+03 y=0.000000000e+00 z=0.000000000e+00
  .default
  +  w=7.500000000e+01 h=2.500000000e+01
  .units mm
  N_5 x=1.500000000e+04 y=0.000000000e+00 z=0.000000000e+00
  G_6
  + x1=1.100000000e+04 y1=1.000000000e+03 z1=-1.000000000e+03
  + x2=1.100000000e+04 y2=9.000000000e+03 z2=-1.000000000e+03
  + x3=1.900000000e+04 y3=9.000000000e+03 z3=-1.000000000e+03
  + thick=2.000000000e+02 seg1=20 seg2=20
  + N_7 (1.200000000e+04,2.000000000e+03,-1.000000000e+03)
  + N_8 (1.300000000e+04,3.000000000e+03,-1.000000000e+03)
  + hole point (1.400000000e+04, 4.000000000e+03, -1.000000000e+03)
  + hole rect (1.500000000e+04, 5.000000000e+03, -1.000000000e+03, 1.600000000e+04, 6.000000000e+03, -1.000000000e+03)
  + hole circle (1.700000000e+04, 7.000000000e+03, -1.000000000e+03, 5.000000000e+02)
  .units um
  N_9 x=1.600000000e+07 y=0.000000000e+00 z=0.000000000e+00
  .units in
  N_10 x=6.692913386e+02 y=0.000000000e+00 z=0.000000000e+00
  .units mils
  N_11 x=7.086614173e+05 y=0.000000000e+00 z=0.000000000e+00
  E_12 N_1 N_2
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_13 N_2 N_3
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_14 N_3 N_4
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_15 N_4 N_5
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_16 N_5 N_9
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_17 N_9 N_10
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_18 N_10 N_11
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  """
  testelement(g2,verified)
  g3 = transform(g2,rx(π/4)*ry(π/4)*rz(π/4))
  verified =
  """
  N_1 x=5.500000000e+00 y=-1.610912703e+00 z=9.389087297e+00
  .units m
  N_2 x=6.000000000e+00 y=-1.757359313e+00 z=1.024264069e+01
  .units km
  N_3 x=6.500000000e-03 y=-1.903805922e-03 z=1.109619408e-02
  .units cm
  N_4 x=7.000000000e+02 y=-2.050252532e+02 z=1.194974747e+03
  .default
  +  w=7.500000000e+01 h=2.500000000e+01
  .units mm
  N_5 x=7.500000000e+03 y=-2.196699141e+03 z=1.280330086e+04
  G_6
  + x1=6.707106781e+03 y1=-1.257359313e+03 z1=8.742640687e+03
  + x2=1.070710678e+04 y2=5.571067812e+03 z2=7.571067812e+03
  + x3=1.470710678e+04 y3=4.399494937e+03 z3=1.439949494e+04
  + thick=2.000000000e+02 seg1=20 seg2=20
  + N_7 (7.707106781e+03,-5.502525317e+02,9.449747468e+03)
  + N_8 (8.707106781e+03,1.568542495e+02,1.015685425e+04)
  + hole point (9.707106781e+03, 8.639610307e+02, 1.086396103e+04)
  + hole rect (1.070710678e+04, 1.571067812e+03, 1.157106781e+04, 1.170710678e+04, 2.278174593e+03, 1.227817459e+04)
  + hole circle (1.270710678e+04, 2.985281374e+03, 1.298528137e+04, 5.000000000e+02)
  .units um
  N_9 x=8.000000000e+06 y=-2.343145751e+06 z=1.365685425e+07
  .units in
  N_10 x=3.346456693e+02 y=-9.801544724e+01 z=5.712758913e+02
  .units mils
  N_11 x=3.543307087e+05 y=-1.037810618e+05 z=6.048803555e+05
  E_12 N_1 N_2
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_13 N_2 N_3
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_14 N_3 N_4
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_15 N_4 N_5
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_16 N_5 N_9
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_17 N_9 N_10
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_18 N_10 N_11
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  """
  testelement(g3,verified)
end
transformtest1()

function transformtest2()
  n0 = Node( 0, 0, 0)
  n1 = Node( 1, 0, 0)
  n2 = Node(-0.1, 0, 0)
  n3 = Node( 0, 0.001, 0)
  n4 = Node( 0,-1.0e-6, 0)
  n5 = Node( 0, 0, 1000)
  n6 = Node( 0, 0,-0.1)
  g4 = Group([
    Units("mm"),
    n0,
    n1,
    Units("cm"),
    n2,
    Units("m"),
    n3,
    Units("km"),
    n4,
    Units("um"),
    n5,
    Units("cm"),
    n6,
    Units("m"),
    Default(w=0.5e-3,h=0.1e-3),
    Segment(n0,n1),
    Segment(n0,n2,wx=0,wy=1,wz=1),
    Segment(n0,n3),
    Segment(n0,n4,wx=1,wy=0,wz=1),
    Segment(n0,n5),
    Segment(n0,n6,wx=1,wy=1,wz=0),
    Units("mm"),
    UniformPlane(x1=2, y1=2, z1=-1.5,
                 x2=2, y2=-2, z2=-1.5,
                 x3=-2, y3=-2, z3=-1.5,
                 thick = .2,
                 seg1 = 20,
                 seg2 = 20)
    ])
  verified =
  """
  .units mm
  N_1 x=0.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  N_2 x=1.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  .units cm
  N_3 x=-1.000000000e-01 y=0.000000000e+00 z=0.000000000e+00
  .units m
  N_4 x=0.000000000e+00 y=1.000000000e-03 z=0.000000000e+00
  .units km
  N_5 x=0.000000000e+00 y=-1.000000000e-06 z=0.000000000e+00
  .units um
  N_6 x=0.000000000e+00 y=0.000000000e+00 z=1.000000000e+03
  .units cm
  N_7 x=0.000000000e+00 y=0.000000000e+00 z=-1.000000000e-01
  .units m
  .default
  +  w=5.000000000e-04 h=1.000000000e-04
  E_8 N_1 N_2
  E_9 N_1 N_3
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=1.000000000e+00
  E_10 N_1 N_4
  E_11 N_1 N_5
  +  wx=1.000000000e+00 wy=0.000000000e+00 wz=1.000000000e+00
  E_12 N_1 N_6
  E_13 N_1 N_7
  +  wx=1.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  .units mm
  G_14
  + x1=2.000000000e+00 y1=2.000000000e+00 z1=-1.500000000e+00
  + x2=2.000000000e+00 y2=-2.000000000e+00 z2=-1.500000000e+00
  + x3=-2.000000000e+00 y3=-2.000000000e+00 z3=-1.500000000e+00
  + thick=2.000000000e-01 seg1=20 seg2=20
  """
  testelement(g4,verified)
  g5 = transform(g4,txyz(10,10,10))
  verified =
  """
  .units mm
  N_1 x=1.000000000e+01 y=1.000000000e+01 z=1.000000000e+01
  N_2 x=1.100000000e+01 y=1.000000000e+01 z=1.000000000e+01
  .units cm
  N_3 x=9.000000000e-01 y=1.000000000e+00 z=1.000000000e+00
  .units m
  N_4 x=1.000000000e-02 y=1.100000000e-02 z=1.000000000e-02
  .units km
  N_5 x=1.000000000e-05 y=9.000000000e-06 z=1.000000000e-05
  .units um
  N_6 x=1.000000000e+04 y=1.000000000e+04 z=1.100000000e+04
  .units cm
  N_7 x=1.000000000e+00 y=1.000000000e+00 z=9.000000000e-01
  .units m
  .default
  +  w=5.000000000e-04 h=1.000000000e-04
  E_8 N_1 N_2
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  E_9 N_1 N_3
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=1.000000000e+00
  E_10 N_1 N_4
  +  wx=-1.000000000e+00 wy=0.000000000e+00 wz=0.000000000e+00
  E_11 N_1 N_5
  +  wx=1.000000000e+00 wy=0.000000000e+00 wz=1.000000000e+00
  E_12 N_1 N_6
  +  wx=1.000000000e+00 wy=0.000000000e+00 wz=0.000000000e+00
  E_13 N_1 N_7
  +  wx=1.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  .units mm
  G_14
  + x1=1.200000000e+01 y1=1.200000000e+01 z1=8.500000000e+00
  + x2=1.200000000e+01 y2=8.000000000e+00 z2=8.500000000e+00
  + x3=8.000000000e+00 y3=8.000000000e+00 z3=8.500000000e+00
  + thick=2.000000000e-01 seg1=20 seg2=20
  """
  testelement(g5,verified)
  g6 = transform(g5,rx(π/4)*ry(π/4)*rz(π/4))
  verified =
  """
  .units mm
  N_1 x=2.928932188e+00 y=1.207106781e+01 z=1.207106781e+01
  N_2 x=3.428932188e+00 y=1.192462120e+01 z=1.292462120e+01
  .units cm
  N_3 x=2.428932188e-01 y=1.221751442e+00 z=1.121751442e+00
  .units m
  N_4 x=3.428932188e-03 y=1.292462120e-02 z=1.192462120e-02
  .units km
  N_5 x=2.428932188e-06 y=1.121751442e-05 z=1.221751442e-05
  .units um
  N_6 x=2.221825407e+03 y=1.257106781e+04 z=1.257106781e+04
  .units cm
  N_7 x=3.636038969e-01 y=1.157106781e+00 z=1.157106781e+00
  .units m
  .default
  +  w=5.000000000e-04 h=1.000000000e-04
  E_8 N_1 N_2
  +  wx=5.000000000e-01 wy=8.535533906e-01 wz=-1.464466094e-01
  E_9 N_1 N_3
  +  wx=-2.071067812e-01 wy=1.353553391e+00 wz=3.535533906e-01
  E_10 N_1 N_4
  +  wx=-5.000000000e-01 wy=1.464466094e-01 wz=-8.535533906e-01
  E_11 N_1 N_5
  +  wx=-2.071067812e-01 wy=3.535533906e-01 wz=1.353553391e+00
  E_12 N_1 N_6
  +  wx=5.000000000e-01 wy=-1.464466094e-01 wz=8.535533906e-01
  E_13 N_1 N_7
  +  wx=1.000000000e+00 wy=7.071067812e-01 wz=7.071067812e-01
  .units mm
  G_14
  + x1=5.989592360e+00 y1=1.273528137e+01 z1=1.273528137e+01
  + x2=3.989592360e+00 y2=9.321067812e+00 z2=1.332106781e+01
  + x3=1.989592360e+00 y3=9.906854249e+00 z3=9.906854249e+00
  + thick=2.000000000e-01 seg1=20 seg2=20
  """
  testelement(g6,verified)
end
transformtest2()


text:
struct Config{T<:Real}
    atol::T # absolute tolerance for ...
    rtol::T # relative tolerance for ...
    solve::Bool # optimize and test result
    query_number_of_constraints::Bool # can get `MOI.NumberOfConstraints` attribute
    query::Bool # can get objective function, and constraint functions, and constraint sets
    modify_lhs::Bool # can modify function of a constraint
    duals::Bool # test dual solutions
    dual_objective_value::Bool # test `DualObjectiveValue`
    infeas_certificates::Bool # check for primal or dual infeasibility certificates when appropriate
    # The expected "optimal" status returned by the solver. Either
    # MOI.OPTIMAL or MOI.LOCALLY_SOLVED.
    optimal_status::MOI.TerminationStatusCode
    basis::Bool # can get variable and constraint basis status

    """
        Config{T}(;
            atol::Real = Base.rtoldefault(T),
            rtol::Real = Base.rtoldefault(T),
            solve::Bool = true,
            query_number_of_constraints::Bool = true,
            query::Bool = true,
            modify_lhs::Bool = true,
            duals::Bool = true,
            dual_objective_value::Bool = duals,
            infeas_certificates::Bool = true,
            optimal_status = MOI.OPTIMAL,
            basis::Bool = false,
        )

    Return an object that is used to configure various tests.

    ## Keywords

     * `atol::Real = Base.rtoldefault(T)`: Control the absolute tolerance used
        when comparing solutions.
     * `rtol::Real = Base.rtoldefault(T)`: Control the relative tolerance used
        when comparing solutions.
     * `solve::Bool = true`: Set to `false` to skip tests requiring a call to
       [`MOI.optimize!`](@ref)
     * `query_number_of_constraints::Bool = true`: Set to `false` to skip tests
       requiring a call to [`MOI.NumberOfConstraints`](@ref).
     * `query::Bool = true`: Set to `false` to skip tests requiring a call to
       [`MOI.get`](@ref) for [`MOI.ConstraintFunction`](@ref) and
       [`MOI.ConstraintSet`](@ref)
     * `modify_lhs::Bool = true`:
     * `duals::Bool = true`: Set to `false` to skip tests querying
       [`MOI.ConstraintDual`](@ref).
     * `dual_objective_value::Bool = duals`: Set to `false` to skip tests
       querying [`MOI.DualObjectiveValue`](@ref).
     * `infeas_certificates::Bool = true`: Set to `false` to skip tests querying
       primal and dual infeasibility certificates.
     * `optimal_status = MOI.OPTIMAL`: Set to `MOI.LOCALLY_SOLVED` if the solver
       cannot prove global optimality.
     * `basis::Bool = false`: Set to `true` if the solver supports
       [`MOI.ConstraintBasisStatus`](@ref) and [`MOI.VariableBasisStatus`](@ref).
    """
    function Config{T}(;
        atol::Real = Base.rtoldefault(T),
        rtol::Real = Base.rtoldefault(T),
        solve::Bool = true,
        query_number_of_constraints::Bool = true,
        query::Bool = true,
        modify_lhs::Bool = true,
        duals::Bool = true,
        dual_objective_value::Bool = duals,
        infeas_certificates::Bool = true,
        optimal_status = MOI.OPTIMAL,
        basis::Bool = false,
    ) where {T<:Real}
        return new(
            atol,
            rtol,
            solve,
            query_number_of_constraints,
            query,
            modify_lhs,
            duals,
            dual_objective_value,
            infeas_certificates,
            optimal_status,
            basis,
        )
    end
    Config(; kwargs...) = Config{Float64}(; kwargs...)
end

@deprecate TestConfig Config

"""
    @moitestset setname subsets

Defines a function `setnametest(model, config, exclude)` that runs the tests defined in the dictionary `setnametests`
with the model `model` and config `config` except the tests whose dictionary key is in `exclude`.
If `subsets` is `true` then each test runs in fact multiple tests hence the `exclude` argument is passed
as it can also contains test to be excluded from these subsets of tests.
"""
macro moitestset(setname, subsets = false)
    testname = Symbol(string(setname) * "test")
    testdict = Symbol(string(testname) * "s")
    if subsets
        runtest = :(f(model, config, exclude))
    else
        runtest = :(f(model, config))
    end
    return esc(
        :(
            function $testname(
                model::$MOI.ModelLike,
                config::$MOI.DeprecatedTest.Config,
                exclude::Vector{String} = String[],
            )
                for (name, f) in $testdict
                    if name in exclude
                        continue
                    end
                    @testset "$name" begin
                        $runtest
                    end
                end
            end
        ),
    )
end


text:
using Statistics

@testset "std" begin
    @test testf(std, rand(10))
    @test testf(std, rand(10,1,2))
    @test testf(std, rand(10,1,2); corrected=true)
    @test testf(std, rand(10,1,2); dims=1)
end

@testset "var" begin
    @test testf(var, rand(10))
    @test testf(var, rand(10,1,2))
    @test testf(var, rand(10,1,2); corrected=true)
    @test testf(var, rand(10,1,2); dims=1)
    @test testf(var, rand(10,1,2); dims=[1])
    @test testf(var, rand(10,1,2); dims=(1,))
    @test testf(var, rand(10,1,2); dims=[2,3])
    @test testf(var, rand(10,1,2); dims=(2,3))
end

@testset "mean" begin
    @test testf(mean, rand(2,2))
    @test testf(mean, rand(2,2); dims=2)
    @test testf(mean, rand(2,2,2); dims=[1,3])
    @test testf(x->mean(sin, x), rand(2,2))
    @test testf(x->mean(sin, x; dims=2), rand(2,2))
    @test testf(x->mean(sin, x; dims=[1,3]), rand(2,2,2))
end


text:
#Demo 10: Mixed-signal "Simulation"
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

line_default = line(color=blue, width=3)
line_overlay = line(color=red, width=3)

#Input data
#-------------------------------------------------------------------------------
ppcycle = 40 #Points per cycle (divisible by 4 to get notches)
T=1/60 #Clock period
ncycles = 15
VSUPPLY = 120 #VRMS - mains
VRECT = 80 #Rectifier voltage
τ = T/2 #Time constant for discharge (output load)

#=="Simulation" (really rudimentary to keep things simple)
===============================================================================#
mains_fn(t) = VSUPPLY*sqrt(2)*cos.(2pi*(t./T))
#Find time where mains crosses given threshold:
mains_xfn(thresh) = acos(thresh/(VSUPPLY*sqrt(2)))*(T/(2pi))

function update_discharge(vrect, Δt, i)
	if 0==vrect[i] #On discharge cycle:
		vrect[i] = VRECT*exp(-Δt/τ)
	end
end

#Compute info for a single cycle:
#-------------------------------------------------------------------------------
t = collect(range(0, stop=T, length=ppcycle+1)[1:end-1]) #Time
npts = length(t); hpts = div(npts, 2) #npts should match ppcycle

mains = mains_fn(t)
mains_rect = abs.(mains)
rectactive = mains_rect .> VRECT
rectactive_sig = rectactive * 1.0
vrect = rectactive_sig*VRECT #Voltage @ rectifier
	#Really dumb algorithm:
	d_start = mains_xfn(VRECT) #Start of 1st discharge cycle
#	d_end = d_start+2*(T/4-d_start)
	for i in 1:hpts
		update_discharge(vrect, t[i]-d_start, i)
	end
	d_start += T/2 #Next discharge cycle
	for i in hpts:npts
		update_discharge(vrect, t[i]-d_start, i)
	end

#Extend data for multiple cycles:
#-------------------------------------------------------------------------------
t = collect(range(0, stop=ncycles*T, length=ncycles*ppcycle+1)[1:end-1]) #Time
mains = repeat(mains, outer=ncycles)
mains_rect = repeat(mains_rect, outer=ncycles)
vrect = repeat(vrect, outer=ncycles)
rectactive_sig = repeat(rectactive_sig, outer=ncycles)

gaterange_list = [4.6:6.1, 9.4:12.9] #In # of cycles
gate_b = 1.0 .+ zero(t) #Active low
for gaterange in gaterange_list
	irange = ppcycle*gaterange
	istart = round(Int, minimum(irange))
	istop = round(Int, maximum(irange))
	gate_b[istart:istop] .= 0
end
vrect_gated = vrect.*gate_b

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Mixed-Signal \"Simulation\"")
mplot.layout[:ncolumns] = 1
ext_digital = InspectDR.PExtents1D(-0.2, 1.2) #"Digital" signals
_E = InspectDR.PExtents1D #Succinct alias for extents.

#To be plotted:
siginfolist = [
	("rectifier", vrect, _E(0,1.2*VRECT)),
	("|mains|>thresh", rectactive_sig, ext_digital),
	("|mains|", mains_rect, _E()),
	("mains", mains, _E()),
]

plot = add(mplot, InspectDR.transientplot([:lin for i in siginfolist],
	title="$(VRECT)V Rectifier",
	ylabels=["(V)" for i in siginfolist])
)
plot.layout[:enable_legend] = true
plot.layout[:halloc_legend] = 150

let wfrm #HIDEWARN_0.7
	for (i, siginfo) in enumerate(siginfolist)
		id, sig, ext = siginfo
		plot.strips[i].yext_full = ext
		wfrm = add(plot, t, sig, id=id, strip=i)
			wfrm.line = line_default
#			wfrm.glyph = InspectDR.glyph(shape=:o, size=3)
	end
end

	#Overlay threshold onto |mains|, for convenience
	wfrm = add(plot, t, zero(t) .+ VRECT, id="thresh", strip=3)
		wfrm.line = line_overlay

	#Overlay gate onto gated output:
	wfrm = add(plot, t, gate_b, id="gate_b", strip=2)
		wfrm.line = line_overlay

	#Overlay gated output onto internal:
	wfrm = add(plot, t, vrect_gated, id="output (gated)", strip=1)
		wfrm.line = line_overlay

gplot = display(InspectDR.GtkDisplay(), mplot)

mplot.layout[:halloc_plot] = 700; mplot.layout[:valloc_plot] = 600
InspectDR.write_png("export_rectifier.png", mplot)

:DONE


text:
export chain_map
"""
    @chain_map(e, v = \_)

The `@chain_map` macro combines three different macros. [`with`](@ref) annotates
each symbol with the associative: `v`. [`chain`](@ref) chains together
expressions wrapped in a `begin` block. [`unweave`](@ref), together with
`NullableArrays.broadcast(lift = true)`, does broadcasting and automatic
lifting of woven expressions.

# Examples
```julia
na =  NullableArrays.NullableArray
a = na(["one", "two"], [false, true])
result = @chain begin
    Dict(:b => na([1, 2]), :c => na(["I", "II"]))
    @chain_map begin
        :b
        sum
        get
        string
        *(~a, " ", _, " ", ~:c)
    end
end

@test get(result[1]) == "one 3 I"
@test result.isnull == [false, true]
```
"""
chain_map(e, v = :_) = @chain begin
    e
    chain
    with(_, v)
    unweave(:(NullableArrays.broadcast(lift = true)), _)
end

@nonstandard chain_map
export @chain_map


text:
function runModels( root::String, modelsDict::Dict, dfd::DataFrame=NA)
   res = OrderedDict()
   q = Symbol[]
   cnt=0
   for (k,v) in OrderedDict(m=>modelsDict[m][:raneff] for m in [:ipen,:idolocc,:iocc])
        for r in v push!(q,k); push!(q,r); cnt=cnt+1 end
   end
   x=reshape(q, (2,cnt))
   ml = permutedims(x, [2, 1])
   m=:empty
   for i in 1:length(ml[:,1])
       if m != ml[i,:][1]
           gr=Symbol("glm_"*string(ml[i,:][1]))
           isna(dfd) ? ( @swarm gr runGlm(root, ml[i,:][1]) ) : ( res[gr] = runGlm(dfd,modelsDict[ml[i,:][1]]) )
           m=ml[i,:][1]
       end
       r=Symbol("glmm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
       isna(dfd) ? (@swarm r runGlmm(root, ml[i,:][1], ml[i,:][2]) ) : (res[r] = runGlmm(dfd,modelsDict[m], [ml[i,:][2]], true) )
       covg=Symbol("covglm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2])   )
       isna(dfd) ? ( @swarm covg runGlm(root, ml[i,:][1], ml[i,:][2]) ) : ( res[covg] = runGlm(dfd,modelsDict[m], ml[i,:][2])  )
   end
   if isna(dfd)
       sleep(20)
       while !iscompletew() println("No!"); sleep(5); end
       println("DONE")
   end
   return res
end

function runModels(root::String, modelsDict::Dict, )
    dfd = read_dfd(root) #dfd = readtable(mod_fname,header=true);
    poolit!(dfd,modelsDict[:factors])
    res = OrderedDict()

    q = Symbol[]
    cnt=0
    for (k,v) in OrderedDict(m=>modelsDict[m][:raneff] for m in [:ipen,:idolocc,:iocc])
         for r in v push!(q,k); push!(q,r); cnt=cnt+1 end
    end
    x=reshape(q, (2,cnt))
    ml = permutedims(x, [2, 1])
    m=:empty
    for i in 1:length(ml[:,1])
        if m != ml[i,:][1]
            gr=Symbol("glm_"*string(ml[i,:][1]))
            res[gr] = runGlm(dfd,modelsDict[ml[i,:][1]]) #res[gr] = runGlm(mod_fname, jmod_fname, ml[i,:][1])
            m=ml[i,:][1]
        end
        r=Symbol("glmm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
        res[r] = runGlmm(dfd,modelsDict[m], [ml[i,:][2]], true)   #res[r] = runGlmm(mod_fname, jmod_fname, ml[i,:][1], ml[i,:][2])
        covg=Symbol("covglm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
        res[covg] = runGlm(dfd,modelsDict[m], ml[i,:][2])   #res[covg] = runGlm(mod_fname, jmod_fname, ml[i,:][1], ml[i,:][2])
    end
    return res
end

function genDHHMeans(dfx::DataFrame)
    println("TOTAL DOLHH")
    o=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="occ")&(dfx[:parameter].=="group"),:adj_mean_score0][1]
    y=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="dolocc")&(dfx[:parameter].=="group"),:adj_mean_score0][1]
    p=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="pen")&(dfx[:parameter].=="group"),:adj_mean_score0][1]
    adj_mean_score0=o*y*p
    println(".... 0 : $adj_mean_score0 ($o, $y,$p) ...")
    o=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="occ")&(dfx[:parameter].=="group"),:adj_mean_score1][1]
    y=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="dolocc")&(dfx[:parameter].=="group"),:adj_mean_score1][1]
    p=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="pen")&(dfx[:parameter].=="group"),:adj_mean_score1][1]
    adj_mean_score1=o*y*p
    println(".... 1 : $adj_mean_score1 ($o, $y,$p) ...")
    unadj_avg_cntrl_hh_pre = mean(dfd[ (dfd[:group] .== 0), :prd_1_net_pr_pre] )
    unadj_avg_expsd_hh_pre = mean(dfd[ (dfd[:group] .== 1), :prd_1_net_pr_pre] )
    unadj_avg_cntrl_hh_pst = mean(dfd[ (dfd[:group] .== 0), :prd_1_net_pr_pos] )
    unadj_avg_expsd_hh_pst = mean(dfd[ (dfd[:group] .== 1), :prd_1_net_pr_pos] )
    pre=["group",0.0,0.0,0.0,0.0,"dolhh","ranef","GLM"]
    vars=[0.0,0.0,adj_mean_score0,adj_mean_score1,unadj_avg_expsd_hh_pre,unadj_avg_expsd_hh_pst,unadj_avg_cntrl_hh_pre,unadj_avg_cntrl_hh_pst]
    cis = zeros(10)  # cis and pvals [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
    pvals=[0.0,0.0]
    cnts=convert(Array{Int64},collect(values(getCnts(dfd))))
    push!(dfx, filler(dfx,vcat(pre,vars,cis,cnts)))
    for ranef in unique(dfx[(dfx[:modelType].=="GLMM"),:ranef])
        for level in unique(dfx[(dfx[:modelType].=="GLMM")&(dfx[:ranef].==ranef),:parameter])
#            println(ranef," ~~ ",level)
            o=dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].=="occ")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:adj_mean_score0][1]
            y=dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].=="dolocc")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:adj_mean_score0][1]
            p=dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].=="pen")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:adj_mean_score0][1]
            adj_mean_score0=o*y*p
            println(".... 1 : $adj_mean_score1 ($o, $y,$p) ...")
            o=dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].=="occ")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:adj_mean_score1][1]
            y=dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].=="dolocc")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:adj_mean_score1][1]
            p=dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].=="pen")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:adj_mean_score1][1]
            adj_mean_score1=o*y*p
            println(".... 0 : ",adj_mean_score0," ... 1 :",adj_mean_score1)
            unadj_avg_expsd_hh_pre = mean(dfd[ (dfd[:group].==1) & (dfd[Symbol(ranef)].==level), :prd_1_net_pr_pre] )
            unadj_avg_expsd_hh_pst = mean(dfd[ (dfd[:group].==1) & (dfd[Symbol(ranef)].==level), :prd_1_net_pr_pos] )
            unadj_avg_cntrl_hh_pre = mean(dfd[ (dfd[:group].==0) & (dfd[Symbol(ranef)].==level), :prd_1_net_pr_pre] )
            unadj_avg_cntrl_hh_pst = mean(dfd[ (dfd[:group].==0) & (dfd[Symbol(ranef)].==level), :prd_1_net_pr_pos] )
            exposed=true # temp hack - all exposed
            if exposed
                unadj_avg_cntrl_hh_pre = dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="dolhh") ,:unadj_avg_cntrl_hh_pre][1]
                unadj_avg_cntrl_hh_pst = dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="dolhh") ,:unadj_avg_cntrl_hh_pst][1]
            end
            pre=[level,0.0,0.0,0.0,0.0,"dolhh",ranef,"GLMM"]
            vars=[0.0,0.0,adj_mean_score0,adj_mean_score1,unadj_avg_expsd_hh_pre,unadj_avg_expsd_hh_pst,unadj_avg_cntrl_hh_pre,unadj_avg_cntrl_hh_pst ]
            cis=zeros(10) #[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
            cnts=convert(Array{Int64},collect(values(getCnts(dfd,ranef,level))))
            push!(dfx, filler(dfx,vcat(pre,vars,cis,cnts)))
        end
    end
end
genDHHMeans(dfx)

ZDict = Dict("onetail_80_pct_intrvl" => 0.84 ,"onetail_90_pct_intrvl" => 1.28, "twotail_80_pct_intrvl" => 1.28, "twotail_90_pct_intrvl" => 1.65)

v_ttl=2

function CIs_O_LB(iDict::OrderedDict, zscore::Float64, iAccuracy::Float64=0.000000001)
    M = get(iDict, :M, NA)
    Mt = get(iDict, :Mt, NA)
    Mc = get(iDict, :Mc, NA)
    N = get(iDict, :N, NA)
    Nt = get(iDict, :Mt, NA)
    Nc = get(iDict, :Nc, NA)
    B1 = get(iDict, :B1, NA)
    B2 = get(iDict, :B2, NA)
    B3 = get(iDict, :B3, NA)
    SE1 = get(iDict, :SE1, NA)
    SE2 = get(iDict, :SE2, NA)
    SE3 = get(iDict, :SE3, NA)
    o_SE0 = get(iDict, :o_SE0, 0)
    y_SE0 = get(iDict, :y_SE0, 0)
    p_SE0 = get(iDict, :p_SE0, 0)
SEsq=sqrt(SE1^2+SE2^2+SE3^2+o_SE0^2+y_SE0^2+p_SE0^2)
    o_B0 = get(iDict, :o_B0, 0)
    y_B0 = get(iDict, :y_B0, 0)
    p_B0 = get(iDict, :p_B0, 0)

    o_mean_score0 = get(iDict, :o_mean_score0, NA)
    o_mean_score1 = get(iDict, :o_mean_score1, NA)
    y_mean_score0 = get(iDict, :y_mean_score0, NA)
    y_mean_score1 = get(iDict, :y_mean_score1, NA)
    p_mean_score0 = get(iDict, :p_mean_score0, NA)
    p_mean_score1 =get(iDict, :p_mean_score1, NA)
Bsum=(B1+B2+B3)-(o_B0+y_B0+p_B0)

    ztot = Bsum-(zscore*SEsq)
    ######CONFIDENCE INTERVAL - LB ########
    m=nothing
    m = Model(solver=NLoptSolver(algorithm=:LD_MMA, maxtime=v_ttl))
    @variable(m, Bocc <= (B1-o_B0))
    @variable(m, Bdolocc <= (B2-y_B0))
    @variable(m, Bpen <= (B3-p_B0))
    @NLobjective(m, Min, ((((p_mean_score1*(Nt/N))+(p_mean_score0*exp(Bpen)*(Nc/N)))*
                                  ((o_mean_score1*(Mt/M))+(o_mean_score0*exp(Bocc)*(Mc/M)))*
                                  ((y_mean_score1*(Mt/M))+(y_mean_score0*exp(Bdolocc)*(Mc/M)))
                                  )
                                  -(((p_mean_score1*(Nt/N)*exp(-Bpen))+(p_mean_score0*(Nc/N)))*
                                    ((o_mean_score1*(Mt/M)*exp(-Bocc))+(o_mean_score0*(Mc/M)))*
                                    ((y_mean_score1*(Mt/M)*exp(-Bdolocc))+(y_mean_score0*(Mc/M)))
                                    )
                                  )
                               )
    @constraint(m, (0.000000000<= (((Bocc+Bpen+Bdolocc)-ztot))<= iAccuracy))
    status = solve(m)
    mval=getobjectivevalue(m)
    println(mval)
    mval_out=mval/(o_mean_score0*y_mean_score0*p_mean_score0)
    return status == :Optimal ? mval_out : -Inf
end
CIs_O_LB(od, 0.84, 0.000000001)

onetail_80 FAILED : DataFrameRow (row 112)
modelType                 GLMM
model                     dolhh
ranef                     ad_type
parameter                 Desktop
adj_dod_effct             10.02936692284249
*****onetail_80_pct_intrvl_lb  10.955272420022512
onetail_80_pct_intrvl_ub  20.010743178509298

LB onetail_80 (ad_type~Desktop):=  - 0.000000001, DataStructures.
od = OrderedDict(:M=>108092,:Mt=>10810,:Mc=>106004,:N=>594930,:Nc=>584120,:B1=>-0.000998217,:B2=>0.0583124,:B3=>0.066694,:SE1=>0.0201122,:SE2=>0.0193296,:SE3=>0.0301203,:o_SE0=>0,:y_SE0=>0,:p_SE0=>0,:o_B0=>0,:y_B0=>0,:p_B0=>0,:o_mean_score0=>1.24344,:o_mean_score1=>1.2422,:y_mean_score0=>3.3212,:y_mean_score1=>3.52063,:p_mean_score0=>0.181431,:p_mean_score1=>0.193944,:metakey=>"ad_type~Desktop",:Bsum=>0.124008,:onetail_pval=>0.998738,:twotail_pval=>0.997476,:twotail_90_pct_intrvl_lb=>0.0678382,:twotail_90_pct_intrvl_ub=>0.245724,:twotail_80_pct_intrvl_lb=>0.0867394,:twotail_80_pct_intrvl_ub=>0.224721)
0.08208352266101582
            : Confidence Interval : 0.10955272420022512

mDict=collectModels(dfx,"GLMM","ad_type","Desktop")
calcPValue_Opt(mDict)
dfx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100

dfx[(dfx[:modelType].=="GLMM")&(dfx[:ranef].=="ad_type")&(dfx[:parameter].=="Desktop"),:]

            adjctl = md[:o_mean_score0] * md[:y_mean_score0] * md[:p_mean_score0]
            adjexp = md[:o_mean_score1] * md[:y_mean_score1] * md[:p_mean_score1]
            println(k," ~~ ",adjctl,"~",adjexp) #,md)
            sdf[sdf[:key].==k,:adj_mean_cntrl_grp] = adjctl
            sdf[sdf[:key].==k,:adj_mean_expsd_grp] = adjexp
            sdf[sdf[:key].==k,:adj_dod_effct] = ((adjexp - adjctl) ./ adjctl ) * 100

function collectModels(dfx::DataFrame,modelType::String,ranef::String="",level::String="")
    mDict = OrderedDict()
    if modelType=="GLM"
        o=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="occ")&(dfx[:parameter].=="group"),:]
        y=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="dolocc")&(dfx[:parameter].=="group"),:]
        p=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="pen")&(dfx[:parameter].=="group"),:]
    else
        o=dfx[(dfx[:modelType].==modelType)&(dfx[:model].=="occ")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:]
        y=dfx[(dfx[:modelType].==modelType)&(dfx[:model].=="dolocc")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:]
        p=dfx[(dfx[:modelType].==modelType)&(dfx[:model].=="pen")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:]
    end
    mDict[:M]=p[:M][1]
    mDict[:Mt]=p[:Mt][1]
    mDict[:Mc]=p[:Mc][1]
    mDict[:N]=p[:N][1]
    mDict[:Mt]=p[:Nt][1]
    mDict[:Nc]=p[:Nc][1]
    if modelType=="GLM"
        mDict[:B1]=o[:coef][1]
        mDict[:B2]=y[:coef][1]
        mDict[:B3]=p[:coef][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100
        mDict[:SE1]=o[:stderr][1]
        mDict[:SE2]=y[:stderr][1]
        mDict[:SE3]=p[:stderr][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100S
    else
        mDict[:B1]=o[:adj_coef][1]
        mDict[:B2]=y[:adj_coef][1]
        mDict[:B3]=p[:adj_coef][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100
        mDict[:SE1]=o[:adj_stderr][1]
        mDict[:SE2]=y[:adj_stderr][1]
        mDict[:SE3]=p[:adj_stderr][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100S
    end
    mDict[:o_SE0]=0
    mDict[:y_SE0]=0
    mDict[:p_SE0]=0
    mDict[:o_B0]=0
    mDict[:y_B0]=0
    mDict[:p_B0]=0
    mDict[:o_mean_score0]=o[:adj_mean_score0][1]
    mDict[:o_mean_score1]=o[:adj_mean_score1][1]
    mDict[:y_mean_score0]=y[:adj_mean_score0][1]
    mDict[:y_mean_score1]=y[:adj_mean_score1][1]
    mDict[:p_mean_score0]=p[:adj_mean_score0][1]
    mDict[:p_mean_score1]=p[:adj_mean_score1][1]
    return mDict
end

function ConfidenceIntervals(dfx::DataFrame)
   # mDict=collectModels(dfx,"GLM")
#    println("Running CI for : Total")
#    mDict[:metakey] = "Total Campaign"
#    calcPValue_Opt(mDict)
#    CIs_O(mDict)
#    dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),:onetail_pval] = mDict[:onetail_pval]*100
#    dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),:twotail_pval] = mDict[:twotail_pval]*100
#    for k in keys(ZDict)
#        dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),Symbol(k*"_lb")] = mDict[Symbol(k*"_lb")]*100
#        dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),Symbol(k*"_ub")] = mDict[Symbol(k*"_ub")]*100
#    end
    for ranef in ["ad_type"] #unique(dfx[(dfx[:modelType].=="GLMM")&(dfx[:ranef].=="ad_type")&(dfx[:parameter].=="Desktop"),:ranef])
        for level in ["Desktop"] #unique(dfx[(dfx[:modelType].=="GLMM")&(dfx[:ranef].==ranef),:parameter])
            mDict=collectModels(dfx,"GLMM",ranef,level)
            println("Running CI for : $ranef : $level")
            mDict[:metakey] = ranef*"~"*level
            calcPValue_Opt(mDict)
            CIs_O(mDict)
            dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:onetail_pval] = mDict[:onetail_pval]*100
            dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:twotail_pval] = mDict[:twotail_pval]*100
            for k in keys(ZDict)
                dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),Symbol(k*"_lb")] = mDict[Symbol(k*"_lb")]*100
                dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),Symbol(k*"_ub")] = mDict[Symbol(k*"_ub")]*100
            end
        end
    end
end
ConfidenceIntervals(dfx)

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

#Generate GLM Formula
#function genFmula(iv::Array{Symbol},isLastOffset::Bool=false)  # Y=1st, if isLastOffset:: offset last
#    isLastOffset ?  genFmula(iv[1],iv,last(iv)) : genFmula(iv[1],iv)
#end
#function genFmula(y::Symbol, iv::Array{Symbol},offset::Symbol=:Missing)
#    xvars=setdiff(iv,[y,offset])
#    eval(parse( string(y)*" ~ 1"* reduce(*, [ " + "*  string(c) for c in xvars ] )  ))
#end
#function genFmula(y::Symbol, iv::Array{Symbol})
#    xvars=setdiff(iv,[y])
#    eval(parse( string(y)*" ~ 1"* reduce(*, [ " + "*  string(c) for c in xvars ] )  ))
#end
function genF(y::Symbol, iv::Array{Symbol},ranef::Array{Symbol}=Symbol[])
        vars=setdiff(iv,vcat([y],ranef))
        f = string(y)*" ~ 1"
        f = length(vars) > 0 ?  f * reduce(*, [ " + "*  string(c) for c in unique(vars) ] ) : f
        f = length(ranef) > 0 ?  f * reduce(*, [ " + "*  "(1 | "*string(c)*")" for c in unique(ranef) ] )  : f
        eval(parse(f))
end

# remove single level vars
function FS_singleLevel(dfd::DataFrame,vars::Array{Symbol}=Symbol[])
    if length(vars)==0  vars=names(dfd)  end
    return [c for c in filter(x -> length(unique( dfd[x] ))==1, vars)]
end

function raneffect(m::MixedModels.GeneralizedLinearMixedModel)
    DataFrame(parameter=levels(m.LMM.trms[1].f), coef=vec(ranef(m, named=true)'[1]), stderr = sqrt.(vec(condVar(m)[1])) )  # stderrX =vec(condVar(m)[1])
end

isnumeric{T<:Number}(::AbstractArray{T}) = true
isnumeric(::Any) = false
isfloat{T<:Float64}(::AbstractArray{T}) = true
isfloat(::Any) = false
isint{T<:Int64}(::AbstractArray{T}) = true
isint(::Any) = false
isString{T<:String}(::AbstractArray{T}) = true
isString(::Any) = false
isBool{T<:Bool}(::AbstractArray{T}) = true
isBool(::Any) = false
isSymbol{T<:Symbol}(::AbstractArray{T}) = true
isSymbol(::Any) = false

Base.lowercase(df::DataFrame) = names!(df, map(x->Symbol(lowercase(string(x))),names(df)))
Base.lowercase(da::DataArray{String}) = [lowercase(x) for x in da]

function filler(dfx::DataFrame, rowIn::Array{Any})
    row=filler(dfx)
    for i in 1:length(rowIn)
        row[i] = rowIn[i]
    end
    return row
end

function filler(dfx::DataFrame, cols::Vector{Symbol}=Symbol[])
    a = Any[]
    if length(cols)==0 cols=names(dfx) end
    for (n, v) in eachcol(dfx)
        if n ∈ cols && isfloat(v)
            push!(a, 0.0)
        elseif n ∈ cols && isint(v)
            push!(a, 0)
        elseif n ∈ cols && isString(v)
            push!(a, "")
        elseif n ∈ cols && isBool(v)
            push!(a, false)
        end
    end
    return a
end

function cov2cor(m)
    d = Diagonal([inv(sqrt(m[i, i])) for i in 1:size(m, 2)])
    A_mul_B!(d, m * d)
end

function vif2(m)
    v = vcov(m)
    assign = m.mm.assign
    nms = coefnames(m.mf)
    if (ipos = findfirst(nms, "(Intercept)")) > 0
        inds = deleteat!(collect(eachindex(nms)), ipos)
        v = view(v, inds, inds)
        assign = view(assign, inds)
    else
        warning("No intercept: vifs may not be sensible.")
    end
    terms = convert(Vector{Symbol}, m.mf.terms.terms)
    n_terms = length(terms)
    if n_terms < 2
        error("model contains fewer than 2 terms")
    end
    R = cov2cor(v)
    detR = det(Symmetric(R))
    result = NamedArray(Array(eltype(R), (length(terms), 3)))
    setnames!(result, string.(terms), 1)
    setnames!(result, ["GVIF", "Df", "GVIF^(1/(2*Df))"], 2)
    inds = collect(eachindex(terms))
    ainds = collect(eachindex(assign))
    for i in inds
        subs = find(x -> x == i, assign)
        notsubs = setdiff(ainds, subs)
        result[i, 1] = det(Symmetric(R[notsubs, notsubs])) / detR
        if (result[i, 2] = length(subs)) > 1
            result[i, 1] *= det(Symmetric(R[subs, subs]))
        end
    end
    if all(result[:, 2] .== 1)
        return result[:, 1:1]
    end
    result
end

function vifDF2(m::RegressionModel)
    sdf=coefDF(m,true)
    v=vif2(m)
    vdf=DataFrame(parameter=convert(Array{Symbol},names(v)[1]), vif=convert(Array{Float64}, v)[:,1]   )
    #vdf[:parameter] = convert(Array{Symbol}, vdf[:variable])
    join(sdf,vdf[[:parameter,:vif]], on = :parameter, kind=:outer)
end

""" Doug
    profileθ(m::GeneralizedLinearMixedModel)

Profile the deviance as a function of θ for models with a single, scalar r.e. term
"""
function profileθ{T}(m::GeneralizedLinearMixedModel{T})
    θ = getθ(m)
    isa(θ, Vector{T}) && length(θ) == 1 ||
        throw(ArgumentError("m must have a single scalar r.e. term"))
    th = zeros(T, 1)
    δ = zeros(T, 1)
    fill!(θ, 0)
    decreasing = false
    basedev = pirls!(setθ!(m, θ), true)
    for t in 0.001:0.001:2.0
        θ[1] = t
        dd = pirls!(setθ!(m, θ), true) - basedev
        sgn = sign(dd - δ[end])
        push!(δ, dd)
        push!(th, t)
        dd ≤ 4.0 || break
        decreasing && sgn > 0 && break
        decreasing = sgn < 0
    end
    th, δ
end

function gethostworkers()
    d = Dict()
    @sync @async for (idx, pid) in enumerate(workers())
        d[pid] = remotecall_fetch(getipaddr,pid)
    end
    iout=Dict()
    for ip in [string(ip) for ip in unique(values(d))]
        a=Int64[]
        for (key, value) in d
            if string(value) == ip
                push!(a,key)
            end
        end
        iout[ip] = sort(a)
    end
    l = get(iout,string(getipaddr()),[])
    iout[string(getipaddr())] = convert(Array{Int64},vcat([1],l))
    return iout
end

function gethostModelDistribution(raneff::Array{Symbol})
    function allocateCore(a::Array{Int64})
        proc = a[1]
        ao = length(setdiff(a,proc)) > 0 ?  setdiff(a,proc) : a
        return proc, ao
    end
    nodehosts = gethostworkers()
    hosts = [k for k in keys(nodehosts)]
    hostcnt = length(hosts)
    #raneff = [:creative, :publisher1, :placement]
    h=1
    m=:pen
    penprocs = Dict()
    for r in raneff
        proc, nodehosts[hosts[h]] = allocateCore(nodehosts[hosts[h]])
        penprocs[r] = Dict(:proc=>proc, :host=>hosts[h])
        println(r,"  :  ", hosts[h], "  ",proc)
        h = h == hostcnt ? h=1 : h+1
    end

    static_nodehosts = deepcopy(nodehosts)

    println("Pen : ",penprocs )
    println(nodehosts)
end

function checksingularity(form::Formula, data::DataFrame, tolerance = 1.e-8)
    mf = ModelFrame(form, data)
    mm = ModelMatrix(mf)
    qrf = qrfact!(mm.m, Val{true})
    vals = abs.(diag(qrf[:R]))
    firstbad = findfirst(x -> x < min(tolerance, 0.5) * vals[1], vals)
    if firstbad == 0
        return Symbol[]
    end
    mf.terms.terms[view(mm.assign[qrf[:p]], firstbad:length(vals))]
end

function coefDF(m::RegressionModel, paramsAsSymbols::Bool=false)
    vout = DataFrame( #vars=vcat([:intercept],m.mf.terms.terms)  #g.model.mm.assign # m.mm.assign
               #vars=vcat([:intercept],[m.mf.terms.terms[x] for x in m.mm.assign[2:end]])
               parameter=DataFrames.coefnames(m.mf)
             , coef=GLM.coef(m)
             , stderr=stderr(m)
             , zval=GLM.coef(m)./stderr(m)
             , pval= ccdf(FDist(1, dof_residual(m)), abs2(GLM.coef(m)./stderr(m))  )
             )
    if paramsAsSymbols == true
        vout[:parameter] = map(x->Symbol(replace(replace(x,"group: 1","group"),"(Intercept)","Intercept")), vout[:parameter])
    end
    return vout
end

"""
Relevel DF for GLM
"""
function relevel{T,I}(v::PooledDataVector{T,I}, lev::T)
    pool = v.pool
    i = convert(I, findfirst(pool, lev))
    i ≠ 0 || throw(ArgumentError("\"$lev\" is not in v.pool"))
    if i == 1
        return copy(v)
    end
    orig = convert(Vector{I}, 1:length(pool))
    perm = append!([i], deleteat!(orig, i))
    PooledDataArray(DataArrays.RefArray(invperm(perm)[v.refs]), v.pool[perm])
end
#dfd[ranef] = relevel(dfd[ranef], "none")

#GLM REsiduals
function xResiduals(g::DataFrames.DataFrameRegressionModel)
    resp = g.model.rr
    #sign(resp.y - resp.mu) .* sqrt(resp.devresid)    # Original
    #sign(resp.y - resp.mu) .* sqrt(complex(resp.devresid))
    sign(resp.y - resp.mu) .* sqrt(((resp.devresid)+1)-1)
    #sign(resp.y - resp.mu) .* sqrt(abs(resp.devresid))
end

function unpool(dfd::DataFrame)
    v_factors = Symbol[]
    for c in names(dfd)
        if typeof(dfd[c]) in [ PooledDataArray{String,UInt8,1} ]
            println("         Converting : ", c )
            push!(v_factors,c)
            dfd[c] = Array(dfd[c])
        end
    end
    return v_factors
end

function poolit!(dfd::DataFrame, vars::Array{Symbol}=Symbol[])
    if length(vars) == 0   vars=names(dfd) end
    for (i, c) in enumerate(dfd.columns)
        if typeof(c) == DataVector{String}
            if dfd.colindex.names[i] in vars
                #println("[",i,"]  : ",dfd.colindex.names[i])
                #unique(dfd[i])
                dfd[i] = pool(dfd[i])
            end
        end
        #typeof(c) == DataVector{String} && println("[",i,"]  : ",dfd.colindex.names[i])
    end
end

function getColswithType(dtype::String, dfd::DataFrame, vars::Array{Symbol}=Symbol[])
    clst=Symbol[]
    typLst = Dict("num" => [DataVector{Float64},DataVector{Int64}], "str" => [DataVector{String}] )
    if length(vars) == 0   vars=names(dfd) end
    if dtype in ["num","str"]
        for (i, c) in enumerate(dfd.columns)
            if dfd.colindex.names[i] in vars
                    if typeof(c) in typLst[dtype]
                #if typeof(c) in [ DataVector{Float64}, DataVector{Int64} ] # = = DataVector{String}
                    #println("[",i,"]  : ",dfd.colindex.names[i])
                    push!(clst,dfd.colindex.names[i])
                end
             end
         end
    else
        println("ERROR: getColswithType : Invalid dtype!")
    end
    return clst
end

function corDFD(dfd::DataFrame, vars::Array{Symbol}=Symbol[])
    if length(vars) == 0
        vars=names(dfd)
    #else
    #    vars=convert(Array{Symbol},vars)
    end
    cols=getColswithType("num", dfd, convert(Array{Symbol},vars) )
    carr= cor(Array(dfd[cols]))
    tcarr = triu(carr,1)
    d = DataFrame(tcarr)
    names!(d,cols)
    d[:vars] = cols
    return stack(d,cols)
end
#function corDFD(dfd::DataFrame, vars::Array{Any}=[])
#    varsS=Symbol[]
#    if length(vars) > 0 varsS=convert(Array{Symbol},vars) end
#    corrDF(dfd,varsS)
#end

# ***************************************************************************************************************************
#export JULIA_NUM_THREADS=4
ENV["JULIA_PKGDIR"] = "/mapr/mapr04p/analytics0001/analytic_users/jpkg"
for i in ["10.106.128.212","10.106.128.213","10.106.128.214","10.106.128.215","10.106.128.216","10.106.128.217",
          "10.106.128.218","10.106.128.219","10.106.128.220","10.106.128.221"]
             if i != string(getipaddr()) println(i); addprocs([(i, 3)]) end
end #addprocs(1)
@everywhere insert!(Base.LOAD_CACHE_PATH, 1, "/mapr/mapr04p/analytics0001/analytic_users/jpkg/lib/v0.5")
@everywhere pop!(Base.LOAD_CACHE_PATH)
using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib, NamedArrays, JuMP, NLopt, NLsolve
@everywhere using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib, NamedArrays, JuMP, NLopt, NLsolve

root = !isdefined(:root) ? pwd() : root
# ************ CUSTOMIZE *********************
#root="/mnt/resource/analytics/models/Jenny-o"
#root="/mnt/resource/analytics/models/NatChoice"
#root="/mnt/resource/analytics/models/rev"
#root="/mnt/resource/analytics/models/ALL#10"
#root="/mnt/resource/analytics/models/ALL#11"
#root="/mnt/resource/analytics/models/CDW#6"
#root="/mnt/resource/analytics/models/HormelChili#8"
#root="/mnt/resource/analytics/models/Rev#8"
#root="/mnt/resource/analytics/models/Jenny-o"
#root="/mnt/resource/analytics/models/hunts"
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/rev"
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/NatChoice"
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/ALL10"

# ************ END CUSTOMIZE *********************

#jmod_fname = root*"/dfd_model.json"
#mod_fname = root*"/dfd_model.csv"
modelsDict = read_modelsDict(root) #modelsDict = readModels(jmod_fname)

@everywhere function runGlm(root::String, modelname::Symbol, ranef::Symbol=:empty)
    println("Loding data : ")
    dfd = read_dfd(root) #dfd = readtable(mod_fname,header=true);
    modelsDict = read_modelsDict(root)#modelsDict = readModels(jmod_fname)
    poolit!(dfd,modelsDict[:factors])
    m=modelsDict[modelname]
    return runGlm(dfd, m, ranef)
end

@everywhere function runGlm(dfd::DataFrame,m::Dict, ranef::Symbol=:empty)
    if ranef==:empty
        f=genF(m[:y_var],m[:finalvars])
        cols = convert(Array{Symbol},vcat(m[:finalvars], [m[:y_var], m[:logvar]]))
        println("Ok - Running ",m[:modelName]," -  Glm on proc ",myid()," with : ",f)
    else
        f=genF(m[:y_var],setdiff(vcat(m[:finalvars],ranef),[:group]))
        dfd[ranef] = relevel(dfd[ranef], "none")
        cols = convert(Array{Symbol},vcat(m[:finalvars], [m[:y_var], m[:logvar]], ranef ))
        println("Ok - Running ",m[:modelName]," - ",ranef," Glm on proc ",myid()," with : ",f)
    end
    if m[:Buyer_Pos_P1_is1]
        return glm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk])
    else
        return return glm(f, dfd[cols] , m[:dist], m[:lnk])
    end
end

@everywhere function runGlmm(root::String, modelname::Symbol, ranef::Symbol=:empty, runθ::Bool=true)
    println("Loding data : Glmm : ",modelname," ~~~ ",ranef)
    dfd = read_dfd(root)   #dfd = readtable(mod_fname,header=true);
    modelsDict = read_modelsDict(root) #modelsDict = readModels(jmod_fname)
    poolit!(dfd,modelsDict[:factors])
    m=modelsDict[modelname]
    if ranef==:empty
        v_out=Dict()
        for r in m[:raneff]
            v_out[r] = runGlmm(dfd, m, [r], runθ)
        end
        return v_out
    else
        return runGlmm(dfd, m, [ranef], runθ)
    end
end

@everywhere function runGlmm(dfd::DataFrame,m::Dict, ranef::Array{Symbol}=[:empty], runθ::Bool=true)
    ranef = ranef==[:empty] ? m[:raneff]  : ranef
    for r in ranef
        dfd[r] = relevel(dfd[r], "none")
    end
    f=genF(m[:y_var],setdiff(m[:finalvars],[:group]),ranef)
    cols = convert(Array{Symbol},vcat(m[:finalvars], [m[:y_var], m[:logvar]], ranef ))
    println("Ok - Running ",m[:modelName]," - ",ranef," Glmm on proc ",myid()," with : ",f)
    if m[:Buyer_Pos_P1_is1]
        #gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), false, 0)   #gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), nAGQ=0)
        gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), verbose=false, fast=true)
        if runθ if sum(raneffect(gmm1)[:coef]) == 0 profileθ(gmm1) end  end
        return gmm1
    else
        #gmm1 = fit!(glmm(f, dfd[cols], m[:dist], m[:lnk]), false, 0)   #gmm1 = fit!(glmm(f, dfd[cols], m[:dist], m[:lnk]), nAGQ=0)
        gmm1 = fit!(glmm(f, dfd[cols] , m[:dist], m[:lnk]), verbose=false, fast=true)
        if runθ if sum(raneffect(gmm1)[:coef]) == 0 profileθ(gmm1) end end
        return gmm1
    end
end

"""
SWARMING - hold on
"""
const SHELF = OrderedDict()
freew() = setdiff(workers(),[v[:worker] for (k,v) in SHELF])
hasfreew() = length(freew()) > 0 ? true : false
nextFreew() = length(freew()) > 0 ? freew()[1]   : NA
#statusW() =  [ string(k)*"   status: "*string(v[:status]) for (k,v) in values(SHELF)]
#statusw() =  for (k,v) in SHELF println(string(k)*"   status: "*string(v[:status])*"   "*string(v[:worker]) ) end
notcompletew() = filter((k,v)-> v[:status]!=:complete ,SHELF)
iscompletew() = length(notcompletew()) == 0 ? true : false
keysw() = collect(keys(SHELF))
getw(k::Symbol) = SHELF[k][:results]
function takew(k::Symbol)
    vo = deepcopy(SHELF[k][:results])
    delete!(SHELF, k)
    return vo
end

type statusW function statusW()  return new() end  end
function Base.show(io::IO, statusw::statusW)  for (k,v) in SHELF println(string(k)*"   status: "*string(v[:status])*"   "*string(v[:worker]) ) end end
statusw=statusW()

function swarm(k::Symbol, thunk)
    #println("SWARMING  : ",k)
    #println("SWARM DUMP : ",typeof(thunk)==Expr)
    SHELF[k] = Dict{Symbol,Any}(:worker =>0,:expr => thunk, :status =>:ready, :channel => Channel(1))
    t = SHELF[k]
    while !hasfreew() sleep(15) end   # println("waiting for free worker");
    t[:worker] = nextFreew()
    wid = t[:worker]
    if isna(wid)
        return false
    else
        #remotecall_fetch(runGlmm, wid, mod_fname, jmod_fname, w[:model], w[:renef] )
        t[:status] = :running
        #splice!(expr.args, 1:1,[:remotecall_fetch,expr.args[1],wid])
        #println("Starting model")
        println("SWARMING  : ",k, " on ", wid)
        @async put!(t[:channel], fetch(Base.sync_add(remotecall(thunk, wid)))  )
        while isready(t[:channel])==false sleep(15) end  #println("NOT Ready : ",t[:worker] );
        t[:results] = take!(t[:channel])
        t[:status] = :complete
        t[:worker] = 0
        println("SWARM Task Completed : ",k)
        close(t[:channel])
        return true
    end
end

macro swarm(p, expr)
    expr = Base.localize_vars(esc(:(()->($expr))), false)
    ex = :( @async swarm($(esc(p)), $expr))
end

# ******************************************

function runModels( root::String, modelsDict::Dict)
   @swarm :glm_ipen runGlm(root, :ipen)
   sleep(25)
   @swarm :glm_iocc runGlm(root, :iocc)
   sleep(25)
   @swarm :glm_idolocc runGlm(root, :idolocc)
   q = Symbol[]
   cnt=0
   for (k,v) in OrderedDict(m=>modelsDict[m][:raneff] for m in [:ipen,:idolocc,:iocc])
        for r in v push!(q,k); push!(q,r); cnt=cnt+1 end
   end
   x=reshape(q, (2,cnt))
   ml = permutedims(x, [2, 1])
   m=:empty
   for i in 1:length(ml[:,1])
       r=Symbol("glmm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
       @swarm r runGlmm(root, ml[i,:][1], ml[i,:][2])
       covg=Symbol("covglm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2])   )
       @swarm covg runGlm(root, ml[i,:][1], ml[i,:][2])
   end
   sleep(20)
   while !iscompletew() println("Not Complete yet!"); sleep(5); end
   println("DONE!DONE!DONE!")
end
runModels(root,modelsDict)

function consolidateResults(modelsDict::Dict,dx::OrderedDict=OrderedDict())
  sdf = DataFrame(parameter=String[], coef=Int64[], stderr=Int64[], zval=Int64[], pval=Int64[], model=String[], ranef=String[], modelType=String[])
    for m in [:iocc, :idolocc, :ipen]
        g = length(dx) > 0 ? dx[Symbol("glm_"*string(m))] : getw( Symbol("glm_"*string(m)) )
        #g = takew(  Symbol("glm_"*string(m))    )
        dfx = coefDF(g,false)
        dfx[:model] = replace(string(m),"i","")
        dfx[:ranef] = "none"
        dfx[:modelType] = "GLM"
        sdf = vcat(sdf,dfx)
        grp = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:coef][1]
        err = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:stderr][1]
        for r in  modelsDict[m][:raneff]
            k="_"*string(m)*"_"*string(r)
            g = length(dx) > 0 ?  dx[Symbol("glmm"*k)]  : getw(Symbol("glmm"*k))
            dfx = raneffect(g)
            dfx[:zval]=0.0
            dfx[:pval]=0.0
            dfx[:model] = replace(string(m),"i","")
            dfx[:ranef]=string(r)
            dfx[:modelType]="GLMM"
            sdf = vcat(sdf,dfx)

            g = length(dx) > 0 ?  dx[Symbol("covglm"*k)]    : getw(Symbol("covglm"*k))
            dfx = coefDF(g,false)
            dfx[:model] = replace(string(m),"i","")
                dfx[:ranef] = string(r)
            dfx[:modelType] = "RanCov"
            dfx = dfx[findin(dfx[:parameter],filter(x->contains(x,string(r)) ,Array(dfx[:parameter]))),:]
            #dfx[:coef]=dfx[:coef]-grp
            #dfx[:stderr]=sqrt(dfx[:stderr].^2+err^2)
            dfx[:parameter] = map(x-> replace(x,string(r)*": ",""),dfx[:parameter])
            sdf = vcat(sdf,dfx)
        end
    end
    return sdf
end
dfx = consolidateResults(modelsDict)
save_dfx(root,dfx)      #writetable(root*"/campaign.csv", dfx)

# ***************************************************************************************************************************
# ***************************************************************************************************************************
# ***************************************************************************************************************************
# ***************************************************************************************************************************
# ***************************************************************************************************************************

#mv(root*"/dfx.csv",root*"/dfx.csv_prev")
#writetable(root*"/dfx_ready.csv", dfx)
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/rev"
#modelsDict = read_modelsDict(root)
#modelsDict[:occ]=modelsDict[:iocc]; delete!(modelsDict,:iocc)
#modelsDict[:dolocc]=modelsDict[:idolocc]; delete!(modelsDict,:idolocc)
#modelsDict[:pen]=modelsDict[:ipen]; delete!(modelsDict,:ipen)
##dfx = read_dfx(root)

#export JULIA_NUM_THREADS=4
ENV["JULIA_PKGDIR"] = "/mapr/mapr04p/analytics0001/analytic_users/jpkg"
for i in ["10.106.128.212","10.106.128.213","10.106.128.214","10.106.128.215","10.106.128.216","10.106.128.217",
          "10.106.128.218","10.106.128.219","10.106.128.220","10.106.128.221"]
             if i != string(getipaddr()) println(i); addprocs([(i, 3)]) end
end #addprocs(1)
@everywhere insert!(Base.LOAD_CACHE_PATH, 1, "/mapr/mapr04p/analytics0001/analytic_users/jpkg/lib/v0.5")
@everywhere pop!(Base.LOAD_CACHE_PATH)
using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib, NamedArrays, JuMP, NLopt, NLsolve
@everywhere using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib, NamedArrays, JuMP, NLopt, NLsolve

"""
SWARMING - hold on
"""
const SHELF = OrderedDict()
freew() = setdiff(workers(),[v[:worker] for (k,v) in SHELF])
hasfreew() = length(freew()) > 0 ? true : false
nextFreew() = length(freew()) > 0 ? freew()[1]   : NA
#statusW() =  [ string(k)*"   status: "*string(v[:status]) for (k,v) in values(SHELF)]
#statusw() =  for (k,v) in SHELF println(string(k)*"   status: "*string(v[:status])*"   "*string(v[:worker]) ) end
notcompletew() = filter((k,v)-> v[:status]!=:complete ,SHELF)
iscompletew() = length(notcompletew()) == 0 ? true : false
keysw() = collect(keys(SHELF))
getw(k::Symbol) = SHELF[k][:results]
function takew(k::Symbol)
    vo = deepcopy(SHELF[k][:results])
    delete!(SHELF, k)
    return vo
end
            # NOTE : create function rmW() # remove value from SHELF
            #                        clearW() #remove all values from SHELF
            #                        waitW()  # while !iscompletew() println("Not Complete yet!"); sleep(5); end

type statusW function statusW()  return new() end  end
function Base.show(io::IO, statusw::statusW)  for (k,v) in SHELF println(string(k)*"   status: "*string(v[:status])*"   "*string(v[:worker]) ) end end
statusw=statusW()

function swarm(k::Symbol, thunk)
    #println("SWARMING  : ",k)
    #println("SWARM DUMP : ",typeof(thunk)==Expr)
    SHELF[k] = Dict{Symbol,Any}(:worker =>0,:expr => thunk, :status =>:ready, :channel => Channel(1))
    t = SHELF[k]
    while !hasfreew() sleep(15) end   # println("waiting for free worker");
    t[:worker] = nextFreew()
    wid = t[:worker]
    if isna(wid)
        return false
    else
        #remotecall_fetch(runGlmm, wid, mod_fname, jmod_fname, w[:model], w[:renef] )
        t[:status] = :running
        #splice!(expr.args, 1:1,[:remotecall_fetch,expr.args[1],wid])
        #println("Starting model")
        println("SWARMING  : ",k)
        @async put!(t[:channel], fetch(Base.sync_add(remotecall(thunk, wid)))  )
        while isready(t[:channel])==false sleep(15) end  #println("NOT Ready : ",t[:worker] );
        t[:results] = take!(t[:channel])
        t[:status] = :complete
        t[:worker] = 0
        println("SWARM Task Completed : ",k)
        close(t[:channel])
        return true
    end
end
macro swarm(p, expr)
    expr = Base.localize_vars(esc(:(()->($expr))), false)
    ex = :( @async swarm($(esc(p)), $expr))
end

@everywhere function MP_collectModels(dfx::DataFrame, modelType::String,ranef::String="",level::String="")
    mDict = OrderedDict()
    if modelType=="GLM"
        o=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="occ")&(dfx[:parameter].=="group"),:]
        y=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="dolocc")&(dfx[:parameter].=="group"),:]
        p=dfx[(dfx[:modelType].=="GLM")&(dfx[:model].=="pen")&(dfx[:parameter].=="group"),:]
    else
        o=dfx[(dfx[:modelType].==modelType)&(dfx[:model].=="occ")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:]
        y=dfx[(dfx[:modelType].==modelType)&(dfx[:model].=="dolocc")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:]
        p=dfx[(dfx[:modelType].==modelType)&(dfx[:model].=="pen")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:]
    end
    mDict[:M]=p[:M][1]
    mDict[:Mt]=p[:Mt][1]
    mDict[:Mc]=p[:Mc][1]
    mDict[:N]=p[:N][1]
    mDict[:Nt]=p[:Nt][1]
    mDict[:Nc]=p[:Nc][1]
    if modelType=="GLM"
        mDict[:B1]=o[:coef][1]
        mDict[:B2]=y[:coef][1]
        mDict[:B3]=p[:coef][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100
        mDict[:SE1]=o[:stderr][1]
        mDict[:SE2]=y[:stderr][1]
        mDict[:SE3]=p[:stderr][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100S
    else
        mDict[:B1]=o[:adj_coef][1]
        mDict[:B2]=y[:adj_coef][1]
        mDict[:B3]=p[:adj_coef][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100
        mDict[:SE1]=o[:adj_stderr][1]
        mDict[:SE2]=y[:adj_stderr][1]
        mDict[:SE3]=p[:adj_stderr][1] #fx[:adj_dod_effct] = ((dfx[:adj_mean_score1] .- dfx[:adj_mean_score0]) ./ dfx[:adj_mean_score0] ) *100S
    end
    mDict[:o_SE0]=0
    mDict[:y_SE0]=0
    mDict[:p_SE0]=0
    mDict[:o_B0]=0
    mDict[:y_B0]=0
    mDict[:p_B0]=0
    mDict[:o_mean_score0]=o[:adj_mean_score0][1]
    mDict[:o_mean_score1]=o[:adj_mean_score1][1]
    mDict[:y_mean_score0]=y[:adj_mean_score0][1]
    mDict[:y_mean_score1]=y[:adj_mean_score1][1]
    mDict[:p_mean_score0]=p[:adj_mean_score0][1]
    mDict[:p_mean_score1]=p[:adj_mean_score1][1]
    println("Running CI for : $ranef : $level")
    mDict[:metakey] = ranef*"~"*level
    calcPValue_Opt(mDict)
    CIs_O(mDict)
    return mDict
end

function MP_ConfidenceIntervals(dfx::DataFrame)
                @swarm :total MP_collectModels(dfx, "GLM")
    for ranef in unique(dfx[(dfx[:modelType].=="GLMM"),:ranef])
        for level in unique(dfx[(dfx[:modelType].=="GLMM")&(dfx[:ranef].==ranef),:parameter])
                k = Symbol(ranef*"~"*level)
                @swarm k  MP_collectModels(dfx, "GLMM",ranef,level)
        end
    end
    sleep(20)
    while !iscompletew() println("Not Complete yet!"); sleep(5); end
    mDict = takew(:total)
    dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),:onetail_pval] = mDict[:onetail_pval]*100
    dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),:twotail_pval] = mDict[:twotail_pval]*100
    dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),:pval] = mDict[:pval]
    for k in keys(ZDict)
        dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),Symbol(k*"_lb")] = mDict[Symbol(k*"_lb")]*100
        dfx[(dfx[:model].=="dolhh")&(dfx[:modelType].=="GLM")&(dfx[:parameter].=="group"),Symbol(k*"_ub")] = mDict[Symbol(k*"_ub")]*100
    end
    for ranef in unique(dfx[(dfx[:modelType].=="GLMM"),:ranef])
        for level in unique(dfx[(dfx[:modelType].=="GLMM")&(dfx[:ranef].==ranef),:parameter])
                mDict = takew( Symbol(ranef*"~"*level) )
                dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),:pval] = mDict[:pval]
                for k in keys(ZDict)
                    dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),Symbol(k*"_lb")] = mDict[Symbol(k*"_lb")]*100
                    dfx[(dfx[:model].=="dolhh")&(dfx[:ranef].==ranef)&(dfx[:parameter].==level),Symbol(k*"_ub")] = mDict[Symbol(k*"_ub")]*100
                end
        end
    end
end

MP_ConfidenceIntervals(dfx)

# ***************************************************************************************************************************
# ***************************************************************************************************************************
# ***************************************************************************************************************************
# ***************************************************************************************************************************
# ***************************************************************************************************************************

# *******************************************
function runModels(root::String, modelsDict::Dict)
    dfd = read_dfd(root) #dfd = readtable(mod_fname,header=true);
    poolit!(dfd,modelsDict[:factors])
    res = OrderedDict()
    q = Symbol[]
    cnt=0
    for (k,v) in OrderedDict(m=>modelsDict[m][:raneff] for m in [:ipen,:idolocc,:iocc])
         for r in v push!(q,k); push!(q,r); cnt=cnt+1 end
    end
    x=reshape(q, (2,cnt))
    ml = permutedims(x, [2, 1])
    m=:empty
    for i in 1:length(ml[:,1])
        if m != ml[i,:][1]
            gr=Symbol("glm_"*string(ml[i,:][1]))
            res[gr] = runGlm(dfd,modelsDict[ml[i,:][1]]) #res[gr] = runGlm(mod_fname, jmod_fname, ml[i,:][1])
            m=ml[i,:][1]
        end
        r=Symbol("glmm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
        res[r] = runGlmm(dfd,modelsDict[m], [ml[i,:][2]], true)   #res[r] = runGlmm(mod_fname, jmod_fname, ml[i,:][1], ml[i,:][2])
        covg=Symbol("covglm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
        res[covg] = runGlm(dfd,modelsDict[m], ml[i,:][2])   #res[covg] = runGlm(mod_fname, jmod_fname, ml[i,:][1], ml[i,:][2])
    end
    return res
end
#dx = runModels(root,modelsDict)

function consolidateResults(modelsDict::Dict,dx::OrderedDict)
  sdf = DataFrame(parameter=String[], coef=Int64[], stderr=Int64[], zval=Int64[], pval=Int64[], model=Symbol[], ranef=String[], modelType=String[])
    for m in [:iocc, :idolocc, :ipen]
        g = dx[Symbol("glm_"*string(m))]
        dfx = coefDF(g,false)
        dfx[:model] = replace(string(m),"i","")
        dfx[:ranef] = "none"
        dfx[:modelType] = "GLM"
        sdf = vcat(sdf,dfx)
        grp = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:coef][1]
        err = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:stderr][1]
        for r in  modelsDict[m][:raneff]
            k="_"*string(m)*"_"*string(r)
            g = dx[Symbol("glmm"*k)]
            dfx = raneffect(g)
            dfx[:zval]=0.0
            dfx[:pval]=0.0
            dfx[:model] = replace(string(m),"i","")
            dfx[:ranef]=string(r)
            dfx[:modelType]="GLMM"
            sdf = vcat(sdf,dfx)
            g = dx[Symbol("covglm"*k)]
            dfx = coefDF(g,false)
            dfx[:model] = replace(string(m),"i","")
                dfx[:ranef] = string(r)
            dfx[:modelType] = "RanCov"
            dfx = dfx[findin(dfx[:parameter],filter(x->contains(x,string(r)) ,Array(dfx[:parameter]))),:]
            dfx[:parameter] = map(x-> replace(x,string(r)*": ",""),dfx[:parameter])
            sdf = vcat(sdf,dfx)
        end
    end
    return sdf
end
#dx = runModels(root,modelsDict)
#dfx = consolidateResults(modelsDict, dx)
#save_dfx(root,dfx)   #writetable(root*"/campaign.csv", dfx)

#dx = runModels(root,modelsDict)
#dfx = consolidateResults(modelsDict, dx)
#save_dfx(root,dfx)   #writetable(root*"/campaign.csv", dfx)

 # ***************************************************************************************************************************

#function distributeDataset()
#    #jmod_fname = root*"/dfd_model.json"
#    #mod_fname
#    #ppp="/mnt/resource/analytics/models/"; run(`ls $ppp`)
#    ips = ["10.63.36.22","10.63.36.23"]
#    for ip in ips
#            cmd=`scp $root/dfd_model.* $ip:$root/`
#            println(cmd)
#            run(cmd)
#        #run(`ssh $ip mkdir $root`)
#        #run(`scp $root/dfd_model.* $ip:$root/`)
#    end
#end
#addprocs([("iriadmin@10.63.36.22", 1), ("iriadmin@10.63.36.23", 1)])
#addprocs(2)

function cmparr(a1::Array{Symbol},a2::Array{Symbol},ignore::Array{Symbol}=Symbol[])
    a1=setdiff(a1,ignore)
    a2=setdiff(a2,ignore)
    a1_only=setdiff(a1,a2)
    a2_only=setdiff(a2,a1)
    l1=length(a1_only)
    l2=length(a2_only)
    println(l1,"--",l2)
    l1>l2 ? DataFrame(a1_only=a1_only,a2_only=vcat(a2_only,fill(:x, l1-l2))) : DataFrame(DataFrame(a1_only=vcat(a1_only,fill(:x, l2-l1)),a2_only=a2_only))
end

function cmparr(a1::Array{Any},a2::Array{Symbol},ignore::Array{Symbol}=Symbol[])
    cmparr(convert(Array{Symbol},a1), a2)
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
# Experiment with different ways to compute Huv:

function Huv_alt(l_max,lam)
Huv = zeros(typeof(lam),l_max+3,l_max+1)
clam = cos(lam); slam = sin(lam)
clam2 = clam*clam; clamn = clam; slamn = slam
for u=0:2:l_max+2
  if u == 0
    Huv[1,1]=  2*lam+pi
    Huv[1,2]= -2*clam
    slamn = slam
    v=2
    while v <= l_max
      Huv[1,v+1]= (-2*clam*slamn+(v-1)*Huv[1,v-1])/(u+v)
      println("lam: ",lam," u: ",u," v: ",v," f1: ",-2*clam*slamn/(u+v)," f2: ",(v-1)/(u+v)*Huv[1,v-1]," Huv: ",Huv[1,v+1])
      slamn *= slam
      v+=1
    end
  else
    slamn = slam
    v=0
    Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
    println("lam: ",lam," u: ",u," v: ",v," f1: ",2*clamn*slamn/(u+v)," f2: ",(u-1)/(u+v)*Huv[u-1,v+1]," Huv: ",Huv[u+1,v+1])
    slamn *= slam
    v=1
    slamn *= slam
    v=1
    Huv[u+1,v+1] = -2*clamn*clam2/(u+1)
#    Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
    println("lam: ",lam," u: ",u," v: ",v," f1: ",2*clamn*slamn/(u+v)," f2: ",(u-1)/(u+v)*Huv[u-1,v+1]," Huv: ",Huv[u+1,v+1])
    slamn = slam
    clamn *= clam2
    v=2
    while v <= l_max
#      Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
      Huv[u+1,v+1]= (-2*clamn*slamn+(v-1)*Huv[u+1,v-1])/(u+v)
      println("lam: ",lam," u: ",u," v: ",v," f1: ",-2*clamn*slamn/(u+v)," f2: ",(v-1)/(u+v)*Huv[u+1,v-1]," Huv: ",Huv[u+1,v+1])
      slamn *= slam
      v+=1
    end
#    clamn *= clam2
  end
end
return Huv
end

function Huv_compute(l_max,lam)
Huv = zeros(typeof(lam),l_max+3,l_max+1)
clam = cos(lam); slam = sin(lam)
clam2 = clam*clam; clamn = clam; slamn = slam
for u=0:2:l_max+2
  if u == 0
    Huv[1,1]=  2*lam+pi
    Huv[1,2]= -2*clam
    slamn = slam
    v=2
    while v <= l_max
      Huv[1,v+1]= (-2*clam*slamn+(v-1)*Huv[1,v-1])/(u+v)
      println("lam: ",lam," u: ",u," v: ",v," f1: ",-2*clam*slamn/(u+v)," f2: ",(v-1)/(u+v)*Huv[1,v-1]," Huv: ",Huv[1,v+1])
      slamn *= slam
      v+=1
    end
  else
    slamn = slam
    v = 0
    while v <= l_max
      Huv[u+1,v+1]= (2*clamn*slamn+(u-1)*Huv[u-1,v+1])/(u+v)
      println("lam: ",lam," u: ",u," v: ",v," f1: ",2*clamn*slamn/(u+v)," f2: ",(u-1)/(u+v)*Huv[u-1,v+1]," Huv: ",Huv[u+1,v+1])
      slamn *= slam
      v+=1
    end
    clamn *= clam2
  end
end
return Huv
end

l_max = 10
epsilon = 1e-3
lam = -pi/2+epsilon
Huv = Huv_compute(l_max,lam)
Huv2= Huv_alt(l_max,lam)
println("Huv diff: ",Huv-Huv2)
lam = pi/2-epsilon
Huv = Huv_compute(l_max,lam)
Huv2= Huv_alt(l_max,lam)
println("Huv diff: ",Huv-Huv2)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark