struct Model{T <: AbstractParam, S <: AbstractFixedParam}
    param::T
    fparam::S
    generate_data::Function
end

function make_models()
    R"source('UM459.R')"
    R"x <- get(load('PSA_MLE.res'))"

    R"x_common <- get_all_models(x, FALSE)"
    R"x_mix <- get_all_models(x, TRUE)"

    @rget x_common
    @rget x_mix

    x_common, x_mix
end

function G1Model(n=200)
    param = Param(n, [-0.15 -0.15 -0.15 -0.15], ## xi
                           [0.75 1.0 0.6 0.15 0.0 0.05], ## beta
                           [-0.01 0.1 -0.01], ## mu
                           [0.6 1.3 0.45], ## S
                           zeros(Float64, 2, 3, 3), ##
                           0.28, ## sigma
                           [1.13], [0.065], [0.9 0.33 1.1])

    C = [1.0 0.41 0.19;
         0.41 1.0 0.49;
         0.19 0.49 1.0]

    param.L[1,:,:] = cholesky!(C).L
            
    Model(param,
          FixedParam(n, 2, 3, 3, 6, 3, 15.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G2Model_common(n=200)
    param = Param(n, transform_zero_sum(vcat([5.86 -1.11 -0.80 -1.45; 0.0 0.0 0.0 0.0])),
                  [0.71 0.19 0.07 1.05 0.32 0.45 0.08 0.19 0.09;
                   0.71 0.19 0.07 1.05 0.32 0.45 0.08 0.19 0.09], ## beta
                  [0.0 0.1 0.0;
                   0.05 0.13 0.14], ## mu
                  [0.52 1.08 0.30; 2.27 4.72 1.30], ## S
                  zeros(Float64, 2, 3, 3), ##
                  0.28, ## sigma
                  ##[0.9, 1.2], ##[log(2)/8, log(2)/5],
                  [1.2, 0.8], ## shape
                  [0.01, 0.03], ## scale
                  [-0.05 1.10 0.71;
                   -0.05 1.10 0.71])

    C1 = [1.0 0.42 0.15;
          0.42 1.0 0.63;
          0.15 0.63 1.0]
    
    param.L[1,:,:] = cholesky(C1).L
    param.L[2,:,:] = cholesky(C1).L

    Model(param,
          FixedParam(n, 2, 3, 3, 9, 3, 14.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G3Model_common(n=200)
    param = Param(n,  [-1.71569    0.396632   0.72337    0.417877;
                       0.943128  -0.229865  -0.555035  -0.249195;
                       0.772563  -0.166767  -0.168335  -0.168682], ## xi
                   [0.680831  0.0179264  0.0295299  0.527697  0.216786  0.198365  0.055459  0.0259366  0.0463118;
                    0.680831  0.0179264  0.0295299  0.527697  0.216786  0.198365  0.055459  0.0259366  0.0463118;
                    0.680831  0.0179264  0.0295299  0.527697  0.216786  0.198365  0.055459  0.0259366  0.0463118], ## beta
                   [0.783827   3.39527   1.3682;
                    0.0834099  1.13367  -0.00579376;
                    -0.0223524  1.36281   0.138263], ## mu
                   [0.41662   2.20853   1.44009;
                    0.489405  0.856195  0.142059;
                    0.677583  1.12779   0.255282], ## S
                  zeros(Float64, 3, 3, 3), ##
                  0.2758719120316683, ## sigma
                  [1.2722676092491088, 0.5853348931121277, 1.8478952942493387], ## shape
                  [0.1324342419290664, 0.0026427656627794963, 0.003291411020296823], ## scale
                  [-0.0437378  1.12338  1.06778;
                   -0.0437378  1.12338  1.06778;
                   -0.0437378  1.12338  1.06778]) ## delta

      L1 = [1.0 0.0 0.0;
            -0.283485   0.958977 0.0;
            -0.0586009  0.688782  0.722597]

     L2 = [1.0 0.0 0.0;
           0.349312  0.937007 0.0;      
           -0.026874  0.690872  0.722478]

    L3 = [1.0 0.0 0.0;
          0.516247  0.85644 0.0;
          0.078312  0.666531  0.741353]
    
    param.L[1,:,:] = L1
    param.L[2,:,:] = L2
    param.L[3,:,:] = L3

    Model(param,
          FixedParam(n, 3, 3, 3, 9, 3, 14.74, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G4Model_common(n=200)
    param = Param(n, [1.05065   -0.272717   -0.597822   -0.30393;
                      0.588372  -0.0421882  -0.0987289  -0.280867;
                      0.342734  -0.141391    0.0141266  -0.0628513;
                      -1.98176    0.456296    0.682425   0.647649], ## xi
                   [0.698363  0.0346776  0.0366413  0.555682  0.327393  0.189299  0.0648854  0.117044  0.0360385;
                    0.698363  0.0346776  0.0366413  0.555682  0.327393  0.189299  0.0648854  0.117044  0.0360385;
                    0.698363  0.0346776  0.0366413  0.555682  0.327393  0.189299  0.0648854  0.117044  0.0360385;
                    0.698363  0.0346776  0.0366413  0.555682  0.327393  0.189299  0.0648854  0.117044  0.0360385], ## beta
                   [0.0847991  1.03647  -0.0311116;
                  -0.162945   1.09494  -0.00677467;
                    0.165787   1.54475   0.202036;  
                    0.732935   3.51918   1.54742], ## mu
                  [0.400097  0.775207  0.0983229;
                   0.784702  1.27895   0.286956;
                   0.472082  0.749551  0.155471;
                   0.427116  2.3425    1.4966], ## S
                  zeros(Float64, 4, 3, 3), ##
                  0.2764, ## sigma
                  [0.5481491252817777, 0.5479779849110288, 1.914072045693419, 1.4028548110541126], ## shape
                  [0.0046433585808854615, 0.003836702125103769, 0.008567616582312225, 0.17172971717241955], ## scale
                  [-0.0335354  1.41136  0.766381;
                  -0.0335354  1.41136  0.766381;
                  -0.0335354  1.41136  0.766381;
                   -0.0335354  1.41136  0.766381]) ## delta

    L1 = [1.0 0.0 0.0;
          0.28681  0.957988 0.0;
          -0.16006  0.754256  0.636772]

    L2 = [1.0 0.0 0.0;
          0.551333   0.834285 0.0;
          0.0561526  0.728513  0.682727]

    L3 = [1.0 0.0 0.0;
          0.25514    0.966904 0.0;
          -0.0490336  0.356662  0.932946]

    L4 = [1.0 0.0 0.0;
          -0.29198    0.956424 0.0;
          -0.0816144  0.662468  0.744631]
    
    param.L[1,:,:] = L1
    param.L[2,:,:] = L2
    param.L[3,:,:] = L3
    param.L[4,:,:] = L4

    Model(param,
          FixedParam(n, 3, 3, 3, 9, 3, 14.74, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G5Model_common(n=200)
    param = Param(n, transform_zero_sum(vcat([3.85913  -0.252416  -0.21723   -1.01946;
                                              3.26331  -0.686945  -0.190965  -1.8576; 
                                              2.74887  -0.664152  -0.982845  -0.358021;
                                              -1.13517   0.853032   0.905174  -0.167422;
                                              0.0 0.0 0.0 0.0])),
                  [0.779325  1.05135    0.341292  0.0664238  -0.00193719   0.0159451;
                   0.171515  0.544798   0.223382  0.0819862  -0.0240127   -0.0398766;
                   0.559038  0.841526   1.7512    0.0646566   0.126851    -0.0518065;
                   0.910531  1.24407    1.00238   0.356472   -0.113962     0.0115417;
                   1.06964   3.18616   -0.977044  1.74197    -2.40273      2.02623], ## beta
                   [0.032456   0.0472063  -0.0119529; 
                    -0.207048   0.201388    0.0234077; 
                    -0.0634234  0.113967    0.00832463;
                    0.0245966  0.169543    0.0534571;
                    0.0694684  0.470818    0.24703], ## mu
                  [0.35 0.95 0.22; 0.41 1.11 0.26; 0.043 0.12 0.03; 0.58 1.56 0.37; 1.47 3.96 0.93], ## S
                  zeros(Float64, 2, 3, 3), ##
                  0.27657, ## sigma
                  [1.68, 0.8, 1.77, 1.57, 2.2], [0.15 0.05 0.36 0.66], [0.53 -1.32 0.63; 0.53 -1.32 0.63; 0.53 -1.32 0.63; 0.53 -1.32 0.63; 0.53 -1.32 0.63])

         L1 = [1.0        0.0       0.0;
           -0.283485   0.918914  0.0;     
           -0.0586009  0.670328  0.679233]

    L2 = [1.0       0.0       0.0;
          0.349312  0.904321  0.0;     
          -0.026874  0.657138  0.669897]

    L3 = [1.0       0.0       0.0;
          0.516247  0.828326  0.0;     
          0.078312  0.63311   0.702438]

    param.L[1,:,:] = L1
    param.L[2,:,:] = L2
    param.L[3,:,:] = L3
    param.L[4,:,:] = L4
    
    Model(param,
          FixedParam(n, 2, 3, 3, 6, 3, 15.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G2Model_mix(n=200)
    param = Param(n, transform_zero_sum(vcat([5.82012 -1.30911 -0.308415 -1.09223; 0.0 0.0 0.0 0.0])),
                           [0.71 1.0 0.53 0.1 0.05 -0.02;
                            0.90 1.27 1.42 0.41 0.39 0.14], ## beta
                           [0.0 0.1 0.0;
                            0.05 0.13 0.14], ## mu
                           [0.5 1.0 0.25; 1.5 3.0 0.75], ## S
                           zeros(Float64, 2, 3, 3), ##
                           0.25, ## sigma
                           [1.59, 1.65], [0.19, 0.52], [0.2 -0.15 0.15; 0.2 -0.15 0.15])

    C1 = [1.0 0.43 0.05;
          0.43 1.0 0.60;
          0.05 0.60 1.0]
    
    param.L[1,:,:] = cholesky(C1).L
    param.L[2,:,:] = cholesky(C1).L

    Model(param,
          FixedParam(n, 2, 3, 3, 6, 3, 15.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G3Model_mix(n=200)
    param = Param(n, transform_zero_sum(vcat([5.58913 -1.24644 -0.346802 -1.04714;
                                              5.01629 -1.89479 0.0111542 -2.12278;
                                              0.0 0.0 0.0 0.0])),
                  [0.772001 1.0691 0.473564 0.0844346 0.0449494 -0.0126505;
                   0.151949 0.273523 0.796271 0.0422973 0.0385719 -0.0586332;
                   0.919505 1.29305 1.46563 0.419051 0.385805 0.129744], ## beta
                  [0.0193879 0.0757485 -0.00638583;
                   -0.210521 0.180814 0.0197556;
                   0.0425674 0.120704 0.143396], ## mu
                  [0.38 0.97 0.26; 0.35 0.89 0.24; 1.12 2.85 0.77], ## S
                  zeros(Float64, 3, 3, 3), ##
                  0.279346, ## sigma
                  [1.6, 1.26, 1.65], [0.2, 0.05, 0.53], [0.184 -0.165 0.127; 0.184 -0.165 0.127; 0.184 -0.165 0.127])

    C1 = [1.0 0.202 -0.059;
          0.202 1.0 0.603;
          -0.059 0.603 1.0]
    
    param.L[1,:,:] = cholesky(C1).L
    param.L[2,:,:] = cholesky(C1).L
    param.L[3,:,:] = cholesky(C1).L

    Model(param,
          FixedParam(n, 3, 3, 3, 6, 3, 15.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G4Model_mix(n=200)
    param = Param(n, transform_zero_sum(vcat([4.17036 -0.314482 -0.36471 -0.947192;
                                              3.12245 -0.669232 -0.122069 -1.90321; 
                                              -1.08051 0.83242 0.907163 -0.150338;
                                              0.0 0.0 0.0 0.0])),
                  [0.75541   1.03187    0.438526  0.0679392   0.00342521   0.0026285;
                   0.158113  0.447918   0.643502  0.3680193   0.520754  -0.0517327;
                   0.910455  1.24207    0.9907    0.54634   -1.25036     0.199782;
                   1.06352   3.14072   -0.872009  1.72911    -2.36841      2.00441], ## beta
                  [0.0174885  0.0740613  -0.0061236;
                   -0.248795   0.161989    0.0174643;
                   0.0432951  0.22873    0.1466806;
                   -0.0283931  0.38582     0.258297], ## mu
                  [0.36 0.94 0.21; 0.41 1.07 0.24; 0.59 1.54 0.35; 1.5 3.92 0.90], ## S
                  zeros(Float64, 4, 3, 3), ##
                  0.27656, ## sigma
                  [sqrt(1.68), sqrt(1.06), sqrt(1.57), sqrt(2.20)],
                  [0.15, 0.05, 0.19, 0.36], [0.55 -1.34 0.64;
                                             0.55 -1.34 0.64;
                                             0.55 -1.34 0.64;
                                             0.55 -1.34 0.64])

    C1 = [ 1.0        0.260983  -0.0176565;
           0.260983   1.0        0.645612; 
           -0.0176565  0.645612   1.0]
    
    param.L[1,:,:] = cholesky(C1).L
    param.L[2,:,:] = cholesky(C1).L
    param.L[3,:,:] = cholesky(C1).L
    param.L[4,:,:] = cholesky(C1).L
    
    Model(param,
          FixedParam(n, 2, 3, 3, 6, 3, 15.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G5Model_mix(n=200)
    param = Param(n, transform_zero_sum(vcat([3.85913  -0.252416  -0.21723   -1.01946;
                                              3.26331  -0.686945  -0.190965  -1.8576; 
                                              2.74887  -0.664152  -0.982845  -0.358021;
                                              -1.13517   0.853032   0.905174  -0.167422;
                                              0.0 0.0 0.0 0.0])),
                  [0.779325  1.05135    0.341292  0.0664238  -0.00193719   0.0159451;
                   0.171515  0.544798   0.223382  0.0819862  -0.0240127   -0.0398766;
                   0.559038  0.841526   1.7512    0.0646566   0.126851    -0.0518065;
                   0.910531  1.24407    1.00238   0.356472   -0.113962     0.0115417;
                   1.06964   3.18616   -0.977044  1.74197    -2.40273      2.02623], ## beta
                   [0.032456   0.0472063  -0.0119529; 
                    -0.207048   0.201388    0.0234077; 
                    -0.0634234  0.113967    0.00832463;
                    0.0245966  0.169543    0.0534571;
                    0.0694684  0.470818    0.24703], ## mu
                  [0.35 0.95 0.22; 0.41 1.11 0.26; 0.043 0.12 0.03; 0.58 1.56 0.37; 1.47 3.96 0.93], ## S
                  zeros(Float64, 2, 3, 3), ##
                  0.27657, ## sigma
                  [1.68, 0.8, 1.77, 1.57, 2.2], [0.15 0.05 0.36 0.66], [0.53 -1.32 0.63; 0.53 -1.32 0.63; 0.53 -1.32 0.63; 0.53 -1.32 0.63; 0.53 -1.32 0.63])

    C1 = [1.0       0.289117  0.010405;
          0.289117  1.0       0.650236;
          0.010405  0.650236  1.0]
        
    param.L[1,:,:] = cholesky(C1).L
    param.L[2,:,:] = cholesky(C1).L
    param.L[3,:,:] = cholesky(C1).L
    param.L[4,:,:] = cholesky(C1).L
    param.L[5,:,:] = cholesky(C1).L
    
    Model(param,
          FixedParam(n, 2, 3, 3, 6, 3, 15.0, reshape(repeat(vcat([0.0, 0.2], 0.5:0.5:14.0), inner=n), n, 30), -1.5, 0.0),
          generate_data)
end

function G2Model_spline(n=200)
    param = ParamSpline(n, [-0.15 -0.15 -0.15 -0.15; 0.15 0.15 0.15 0.15], ## xi
                        [2 2 2 0.2 0.3 0.4;
                         1 1 1 0.1 0.1 0.1], ## beta
                        [0 0.0 0.0;
                         0 1.5 0.95], ## mu
                        sqrt.([0.5 1.3 0.3; 0.1 0.1 0.5]), ## S
                        zeros(Float64, 2, 3, 3), ##
                        0.25, ## sigma
                        log.([1.0 1.2 1.1; 0.8 0.6 0.9]))

    param.L[1,:,:] = [1.0 0.0 0.0;
                      0.5 0.866025 0.0;
                      -0.75 0.583124 0.312197]

    param.L[2,:,:] = [1.0 0.0 0.0;
                      0.5 0.866025 0.0;
                      0.75 -0.282902 0.597885]

    param.spline[1] = SplineHazard.Spline.const_spline(log(2)/8, (0.0, 24.0))
    param.spline[2] = SplineHazard.Spline.const_spline(log(2)/5, (0.0, 24.0))
    
    Model(param,
          FixedParam(n, 2, 3, 3, 6, 3, 10.0, reshape(repeat([0.0, 0.5, 1.0, 1.5, 2.0:1.0:10.0...], inner=n), n, 13), -1.5, 0.0),
          generate_data)
end

function rj_mix_test(M::Int; chains=1, warmup=Int(floor(M/2)))

    model2 = JLCMWeibull.G2Model_common(200)
    model3 = JLCMWeibull.G3Model_common(200)
    model4 = JLCMWeibull.G4Model_common(400)
    
    data2 = model2.generate_data(model2.param, model2.fparam)
    data3 = model3.generate_data(model3.param, model3.fparam)
    data4 = model4.generate_data(model4.param, model4.fparam)
    
    data = JLCMWeibull.merge_data(data2, data3)
    
    prior = JLCMWeibull.Prior(3, 3, 6, 3)
    
    x = JLCMWeibull.fit(M, data, model4.param, prior; warmup=warmup, G=3, chains=chains)
    x, data ##, JLCMWeibull.summary(x, data)
end
