text:
@generated function Base.getindex{T}(bag::AbstractMixedBag, ::Type{T})
    :(bag.$(fieldname(T, eltypes(bag))))
end

function Base.push!(bag::AbstractMixedBag, item)
    T = typeof(item)
    push!(bag[T], item)
    bag
end

function Base.append!(bag::AbstractMixedBag, iter)
    T = eltype(iter)
    append!(bag[T], iter)
    bag
end

# comparison
@generated function fieldwise_equality(x, y)
    @assert fieldnames(x) == fieldnames(y)  # better look at typeof(x).types ?
    args = [:(x.$f == y.$f || return false) for f in fieldnames(x)]
    Expr(:block, args..., :(return true))
end

import Base.==
(==)(x::AbstractMixedBag, y::AbstractMixedBag) = fieldwise_equality(x, y)

# collect
function collect_impl(Bag)
    args = map(fname -> :(bag.$(fname)) ,fieldnames(Bag))
    Expr(:call, :vcat, args...)
end
@generated Base.collect(bag::AbstractMixedBag) = collect_impl(bag)

# mapreduce

import Base.mapreduce
@generated mapreduce(f, op, v0, bag::AbstractMixedBag) = mapreduce_impl(f, op, v0, bag)
@generated mapreduce(f, op, bag::AbstractMixedBag) = mapreduce_impl(f, op, bag)

function mapreduce_impl(f, op, v0, Bag)
    fnames = fieldnames(Bag)
    N = length(fnames)
    v(i) = Symbol("v$i")
    args = map((i,fname) -> :($(v(i+1)) = mapreduce(f, op, $(v(i)), bag.$fname) ), 0:(N-1), fnames)
    retarg = :(return $(v(N)))
    Expr(:block, args..., retarg)
end

function mapreduce_impl(f, op, Bag)
    fnames = fieldnames(Bag)
    N = length(fnames)
    v(i) = Symbol("v$i")
    fname1 = fnames[1]
    arg1 = :( $(v(1)) = mapreduce(f, op, bag.$fname1) )
    args = map((i,fname) -> :( $(v(i+1)) = mapreduce(f, op, $(v(i)), bag.$fname) ), 1:(N-1), fnames[2:N])
    retarg = :(return $(v(N)))
    Expr(:block, arg1, args..., retarg)
end

# foreach
Base.foreach{T}(f, bag::AbstractMixedBag, ::Type{T}) = foreach(f, bag[T])
@generated Base.foreach(f, bag::AbstractMixedBag) = foreach_impl(f, bag)
function foreach_impl(f, bag)
    typs = map(eltype, bag.types)
    args = map(T -> :(foreach(f, bag, $T)), typs)
    Expr(:block, args...)
end


text:
@testset "Communications" begin
  server() do port
    @test port isa Int32
    @test hget(("", port), "1 2 3") == [1, 2, 3]
    @test 42 == hopen(port) do h
      hget(h, "42")
    end
    @test_throws KdbException hget(("", port), "(..")
    @test "type" == try
      hget(("", port), "1+`")
    catch e
      e.s
    end
    hopen(port) do h
      auto_r0(k, h, "42") do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{y}", kj(1), kj(42)) do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{z}", kj(1), kj(2), kj(42)) do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{[a;b;c;d]d}", kj(1), kj(2), kj(3), kj(42)) do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{[a;b;c;d;e]e}",
              kj(1), kj(2), kj(3), kj(4), kj(42)) do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{[a;b;c;d;e;f]f}",
              kj(1), kj(2), kj(3), kj(4), kj(5), kj(42)) do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{[a;b;c;d;e;f;g]g}",
              kj(1), kj(2), kj(3), kj(4), kj(5), kj(6), kj(42)) do x
        @test xj(x) == 42
      end
      auto_r0(k, h, "{[a;b;c;d;e;f;g;h]h}",
              kj(1), kj(2), kj(3), kj(4), kj(5), kj(6), kj(7), kj(42)) do x
        @test xj(x) == 42
      end
      @test 0 == hget(h, "0")
      @test 1 == hget(h, "{[a]a}", 1)
      @test 2 == hget(h, "{[a;b]b}", 1, 2)
      @test 3 == hget(h, "{[a;b;c]c}", 1, 2, 3)
      @test 4 == hget(h, "{[a;b;c;d]d}", 1, 2, 3, 4)
      @test 5 == hget(h, "{[a;b;c;d;e]e}", 1, 2, 3, 4, 5)
      @test 6 == hget(h, "{[a;b;c;d;e;f]f}", 1, 2, 3, 4, 5, 6)
      @test 7 == hget(h, "{[a;b;c;d;e;f;g]g}", 1, 2, 3, 4, 5, 6, 7)
      @test 8 == hget(h, "{[a;b;c;d;e;f;g;h]h}", 1, 2, 3, 4, 5, 6, 7, 8)
    end
  end # server()
  server(user="star:shine") do port
    @test 0 < (h = khpu("", port, "star:shine"); hclose(h); h)
    @test 0 < (h = khpun("", port, "star:shine", 10); hclose(h); h)
    @test 42 == hopen(port=port, user="star:shine") do h
      hget(h, "42")
    end
    @test 42 == hopen(port=port, user="star:shine", timeout=1) do h
      hget(h, "42")
    end
  end
end

ext: jl


text:
import Base: !
export !

immutable Unary <: AFImpl
	not::Ptr{Void}

	function Unary(ptr)
		new(
			Libdl.dlsym(ptr, :af_not)
		)
	end
end

macro unOp(op, cFunc)
	:( @afCall_Arr_Arr($(esc(op)), unary, $cFunc) )
end

macro logicUnOp(op, cFunc)
	:( @unOp($(esc(op)), $cFunc) )
end

macro arUnOp(op, cFunc)
	:( @unOp($(esc(op)), $cFunc) )
end

@logicUnOp(!, not)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
println("Alpha-parameter, Cross check 'frohlichalpha()' fn vs. literature values.")

@testset "FrohlichAlpha" begin

α=frohlichalpha(2.3, 5.6, 49 / 2π, 1.0)
println("NaCl Frohlich paper, α=",α," should be ~about 5 (Feynman1955)")
@test α ≈ 5.0 atol=0.3

α=frohlichalpha(7.1, 10.4, 5.08, 0.095)
println("CdTe  α=",α," Stone 0.39 / Devreese 0.29")
@test α ≈ 0.3 atol=0.1

α=frohlichalpha(10.89, 12.9, 8.46, 0.063)
println("GaAs  α=",α," Devreese 0.068 ")
@test α ≈ 0.068 atol=0.01

end


text:
include("Include.jl")

# Evaluates the objective function values -
function local_refienment_step(path_to_data_dir, parameter_array; sigma=0.05, iteration_max=100)

    # initialize -
    number_of_parameters = length(parameter_array)
    BIG = 1e10

    # load the default data_dictionary -
    time_start = 0.0
    time_stop = 16.0
    time_step_size = 0.01

    # what is the host_type?
    host_type = :cell_free

    # path to parameters -
    path_to_biophysical_constants_file = "./CellFree.json"

    # wght array -
    W = diagm(ones(5))
    W[3,3] = 0.1
    W[4,4] = 1000.0
    W[5,5] = 1.0

    # Load the data dictionary (uses the default biophysical_constants file)
    default_data_dictionary = build_data_dictionary((time_start,time_stop,time_step_size), path_to_biophysical_constants_file, host_type)

    # Set some values for the weights, gene lengths etc -
    model_data_dictionary = customize_data_dictionary(default_data_dictionary,host_type)

    # load the experimental data -
    exp_data_dictionary = load_experimental_data_dictionary(path_to_data_dir)

    # setup the functions -
    OF(P) = objective_function(P,time_start,time_step_size,time_stop,model_data_dictionary,exp_data_dictionary)

    function _compute_error_total(objective_array,W)
        value = transpose(objective_array)*W*objective_array
        return value[1]
    end

    # calculate the starting error -
    parameter_array_best = parameter_array
    error_array = BIG*ones(4)
    error_array[1] = _compute_error_total(OF(parameter_array_best), W)

    # main refinement loop -
    iteration_counter = 1
    while (iteration_counter<iteration_max)

        # take a step up -
        parameter_up = parameter_array_best.*(1 .+ sigma*rand(number_of_parameters))
        parameter_up = check_parameter_bounds(parameter_up)

        # take a step down -
        parameter_down = parameter_array_best.*(1 .- sigma*rand(number_of_parameters))
        parameter_down = check_parameter_bounds(parameter_down)

        # Evaluate the obj function -
        error_array[2] = _compute_error_total(OF(parameter_up),W)
        error_array[3] = _compute_error_total(OF(parameter_down),W)

        # Calculate a correction factor -
        a = error_array[2]+error_array[3] - 2.0*error_array[1]
        parameter_corrected = parameter_array_best
        if (a>0.0)
            amda = -0.5*(error_array[3] - error_array[2])/a
            parameter_corrected = parameter_array_best .+ amda*rand(number_of_parameters)
            parameter_corrected = check_parameter_bounds(parameter_corrected)
            error_array[4] = _compute_error_total(OF(parameter_corrected), W)
        end

        # Which step has the min error?
        min_index = argmin(error_array)
        if (min_index == 1)
            parameter_array_best = parameter_array_best
        elseif (min_index == 2)
            parameter_array_best = parameter_up
        elseif (min_index == 3)
            parameter_array_best = parameter_down
        elseif (min_index == 4)
            parameter_array_best = parameter_corrected
        end

        # Update the local error
        error_array[1] = error_array[min_index]

        @show iteration_counter,error_array[min_index]

        # update local counter -
        iteration_counter = iteration_counter + 1
    end

    return parameter_array_best
end

function check_parameter_bounds(parameter_array)

    # setup paramter bounds -
    pvec_bounds = [

        # dG's -
        39995.69565886529   68282.2841495398    ;   # 1 W_clssrA_RNAP
        -50000.0  -1.0   ;   # 2 W_clssrA_sigma_28
        39995.69565886529   68282.2841495398    ;   # 3 W_deGFPssrA_RNAP
        -35689.208748865043  -24993.98186199351    ;   # 4 W_deGFPssrA_sigma_70
        -50000.0  20000.0   ;   # 5 W_deGFPssrA_clssrA
        39995.69565886529   68282.2841495398    ;   # 6 W_sigma_28_RNAP
        -50000.0  25000.0   ;   # 7 W_sigma_28_sigma_70
        -50000.0  25000.0   ;   # 8 W_sigma_28_clssrA

        # binding parameters -
        0.5 10.0            ;   # 9 n_clssrA_sigma_28
        0.001 100.0         ;   # 10 K_clssrA_sigma_28
        0.5037004212911795  2.2       ;   # 11 n_deGFPssrA_sigma_70
        25.539812443045424 100.0         ;   # 12 K_deGFPssrA_sigma_70
        0.5 10.0            ;   # 13 n_deGFPssrA_clssrA
        0.001 100.0         ;   # 14 K_deGFPssrA_clssrA
        0.5 10.0            ;   # 15 n_sigma_28_sigma_70
        0.001 100.0         ;   # 16 K_sigma_28_sigma_70
        0.5 10.0            ;   # 17 n_sigma_28_clssrA
        0.001 100.0         ;   # 18 K_sigma_28_clssrA

        # time constants -
		0.0001 100.0         ;   # 19	mRNA_clssrA
		0.0039273045649573235 7.17667775643900         ;	# 20	mRNA_deGFPssrA
		0.0001 100.0         ;	# 21	mRNA_sigma_28
		0.0001 100.0         ;	# 22	protein_clssrA
		0.35035678021027006 10.203484322453043         ;	# 23	protein_deGFPssrA
		0.001 100.0	        ;	# 24	protein_sigma_28

        # degradation mods -
        0.1 10.0            ;	# 25	mRNA_clssrA 0.0085
		0.0011598053735671 3.6284714410194088 	     ;	# 26	mRNA_deGFPssrA 0.05
		0.1 10.0 	        ;	# 27	mRNA_sigma_28
		5.0 10.0 	        ;	# 28	protein_clssrA
		40.0 60.0 	        ;	# 29	protein_deGFPssrA
		0.13 10.0 	        ;	# 30	protein_sigma_28
        0.132047561240967526 10.137258626014857 	        ;	# 31	protein_sigma_70

        # w -
        4.0 10.0            ;   # 32    translation capacity half-life

        # KL -
        10.81768297340375 396.3396750321431         ;   # 33    KL in muM
    ];

    # tmp -
    pvec_initial = parameter_array

    # check bounds -
    number_of_parameters = length(pvec_initial)
    for parameter_index = 1:number_of_parameters

        # what is the parameter value?
        p_i = pvec_initial[parameter_index]

        # is p_i outside of the bounds?
        lb_value = pvec_bounds[parameter_index,1]
        ub_value = pvec_bounds[parameter_index,2]

        if (p_i<lb_value)
            pvec_initial[parameter_index,1] = lb_value
        end

        if (p_i>ub_value)
            pvec_initial[parameter_index,1] = ub_value
        end
    end

    # return -
    return pvec_initial
end

function neighbor_function(parameter_array; sigma=0.05)

    # setup -
    number_of_parameters = length(parameter_array)

    # calculate new parameter array -
    new_parameter_array = parameter_array.*(1 .+ sigma*randn(number_of_parameters))

    # check the bounds and return -
    return check_parameter_bounds(new_parameter_array)
end

function cooling_function(temperature)

  # define my new temperature -
  alpha = 0.9
  return alpha*temperature
end

function acceptance_probability_function(rank_array,temperature)
    return (exp(-rank_array[end]/temperature))
end

function objective_function(parameter_guess_array,time_start,time_step_size,time_stop,model_data_dictionary,exp_data_dictionary)

    # what is the host_type?
    host_type = :cell_free

    # some global parameters -
    BIG = 1e10
    SMALL = 1e-6

    # Phase 1: parameter update =========================================================================== #
    # update the paramaters in the model data dictionary -
    # for now - lets only search over dG's -
    R = model_data_dictionary["R"]
    T_K = model_data_dictionary["T_K"]

    # compute W -
    tmp_W_array = Float64[]
    for index = 1:8
        parameter_guess = parameter_guess_array[index]
        value = exp(-1*parameter_guess/(R*T_K))
        push!(tmp_W_array,value)
    end

    # update the control W's -
    control_parameter_dictionary = model_data_dictionary["control_parameter_dictionary"]
    control_parameter_dictionary["W_clssrA_RNAP"] = tmp_W_array[1]          # 1
	control_parameter_dictionary["W_clssrA_sigma_28"] = tmp_W_array[2]      # 2
	control_parameter_dictionary["W_deGFPssrA_RNAP"] = tmp_W_array[3]       # 3
	control_parameter_dictionary["W_deGFPssrA_sigma_70"] = tmp_W_array[4]   # 4
	control_parameter_dictionary["W_deGFPssrA_clssrA"] = tmp_W_array[5]     # 5
	control_parameter_dictionary["W_sigma_28_RNAP"] = tmp_W_array[6]        # 6
	control_parameter_dictionary["W_sigma_28_sigma_70"] = tmp_W_array[7]    # 7
	control_parameter_dictionary["W_sigma_28_clssrA"] = tmp_W_array[8]      # 8
    model_data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary

    binding_parameter_dictionary = model_data_dictionary["binding_parameter_dictionary"]
    binding_parameter_dictionary["n_clssrA_sigma_28"] = parameter_guess_array[9]
	binding_parameter_dictionary["K_clssrA_sigma_28"] = parameter_guess_array[10]
	binding_parameter_dictionary["n_deGFPssrA_sigma_70"] = parameter_guess_array[11]
	binding_parameter_dictionary["K_deGFPssrA_sigma_70"] = parameter_guess_array[12]
	binding_parameter_dictionary["n_deGFPssrA_clssrA"] = parameter_guess_array[13]
	binding_parameter_dictionary["K_deGFPssrA_clssrA"] = parameter_guess_array[14]
	binding_parameter_dictionary["n_sigma_28_sigma_70"] = parameter_guess_array[15]
	binding_parameter_dictionary["K_sigma_28_sigma_70"] = parameter_guess_array[16]
	binding_parameter_dictionary["n_sigma_28_clssrA"] = parameter_guess_array[17]
    binding_parameter_dictionary["K_sigma_28_clssrA"] = parameter_guess_array[18]
    model_data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary

    # time constant modifier -
    time_constant_modifier_array = [
		0.0	                          ;	  # 1	clssrA
		0.0	                          ;	  # 2	deGFPssrA
		0.0	                          ;	  # 3	sigma_28
		0.0                           ;	  # 4	sigma_70
		parameter_guess_array[19]     ;	  # 5	mRNA_clssrA
		parameter_guess_array[20]     ;	  # 6	mRNA_deGFPssrA
		parameter_guess_array[21]     ;	  # 7	mRNA_sigma_28
		1.0 	                      ;	  # 8	mRNA_sigma_70
		parameter_guess_array[22]	  ;	  # 9	protein_clssrA
		parameter_guess_array[23]	  ;	  # 10	protein_deGFPssrA
		parameter_guess_array[24]	  ;	  # 11	protein_sigma_28
		1.0	                          ;	  # 12	protein_sigma_70
	]
    model_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

    # setup degradation_modifier_array -
    degradation_modifier_array = [
		0.0	                            ;	# 1	clssrA
		0.0	                            ;	# 2	deGFPssrA
		0.0	                            ;	# 3	sigma_28
		0.0	                            ;	# 4	sigma_70
		parameter_guess_array[25]	    ;	# 5	mRNA_clssrA 0.0085
		parameter_guess_array[26]		;	# 6	mRNA_deGFPssrA 0.05
		parameter_guess_array[27]		;	# 7	mRNA_sigma_28
		0.0                             ;   # 8 mRNA_sigma_70
		parameter_guess_array[28]		;	# 9	protein_clssrA
		parameter_guess_array[29]		;	# 10	protein_deGFPssrA
		parameter_guess_array[30]		;	# 11	protein_sigma_28
	    parameter_guess_array[31]		;	# 12	protein_sigma_70
    ]

    # update the translation time -
    model_data_dictionary["half_life_translation_capacity"] = parameter_guess_array[32]

    # lastly, update KL -
    biophysical_constants_dictionary = model_data_dictionary["biophysical_constants_dictionary"]
    biophysical_constants_dictionary["translation_saturation_constant"] = parameter_guess_array[33]
    model_data_dictionary["biophysical_constants_dictionary"] = biophysical_constants_dictionary

    # grab defaults -
    species_symbol_type_array = model_data_dictionary["species_symbol_type_array"]
    protein_coding_length_array = model_data_dictionary["protein_coding_length_array"]
    gene_coding_length_array = model_data_dictionary["gene_coding_length_array"]
    time_constant_modifier_array = model_data_dictionary["time_constant_modifier_array"]
    initial_condition_array = model_data_dictionary["initial_condition_array"]

    # # get gene IC -
    idx_gene = findall(x->x==:gene,species_symbol_type_array)
    gene_abundance_array = initial_condition_array[idx_gene]

    # Precompute the translation parameters -
	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
    model_data_dictionary["translation_parameter_array"] = translation_parameter_array

	# Precompute the kinetic limit of transcription -
	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
    model_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

    # Dilution degrdation matrix -
    dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)
    model_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix
    # ===================================================================================================== #

    # Phase 2:  solve model equations ===================================================================== #
    # solve the balance equations -
    (TSIM,XSIM) = SolveBalances(time_start,time_stop,time_step_size,model_data_dictionary)
    # ===================================================================================================== #

    # Phase 3:  compute simulation error ================================================================== #
    # compute the error - we need to do a bunch of interpolation -
    tsim_exp = exp_data_dictionary["deGFP_mRNA_data_table"][!,:Time_in_hr]

    # GFP mRNA -
    itp_gfp_mRNA =  LinearInterpolation(TSIM, (1000)*XSIM[:,6]);
    mRNA_GFP_sim = itp_gfp_mRNA[tsim_exp]  # convert to muM from nM

    # sigma 28 mRNA -
    itp_sigma28_mRNA =  LinearInterpolation(TSIM, (1000)*XSIM[:,7]);
    mRNA_sigma28_sim = itp_sigma28_mRNA[tsim_exp] # convert to muM from nM

    # C1 mRNA -
    itp_C1_mRNA =  LinearInterpolation(TSIM, (1000)*XSIM[:,5]);
    mRNA_C1_sim = itp_C1_mRNA[tsim_exp] # convert to muM from nM

    # GFP protein -
    itp_gfp_protein =  LinearInterpolation(TSIM, XSIM[:,10]);
    protein_GFP_sim = itp_gfp_protein[tsim_exp]

    # get experimental data -
    mRNA_GFP_exp = exp_data_dictionary["deGFP_mRNA_data_table"][!,:Concentration_in_nM]
    mRNA_GFP_std_exp = exp_data_dictionary["deGFP_mRNA_data_table"][!,:stdev]

    mRNA_sigma28_exp = exp_data_dictionary["sigma28_mRNA_data_table"][!,:Concentration_in_nM]
    mRNA_sigma28_std_exp = exp_data_dictionary["sigma28_mRNA_data_table"][!,:stdev]

    mRNA_C1_exp = exp_data_dictionary["c1_ssrA_mRNA_data_table"][!,:Concentration_in_nM]
    mRNA_C1_std_exp = exp_data_dictionary["c1_ssrA_mRNA_data_table"][!,:stdev]

    protein_GFP_exp = exp_data_dictionary["deGFP_protein_data_table"][!,:Concentration_in_muM]
    protein_GFP_std_exp = exp_data_dictionary["deGFP_protein_data_table"][!,:stdev]

    # compute error terms -
    error_term_array = BIG*ones(5,1)

    # mRNA GFP -
    #WM1 = diagm(1.0 ./((mRNA_GFP_std_exp).^2))
    error_vector_1 = (mRNA_GFP_exp .- mRNA_GFP_sim)
    error_term_array[1,1] = 1.0*transpose(error_vector_1)*error_vector_1

    # mRNA sigma28 -
    #WM2 = diagm(1.0 ./((mRNA_sigma28_std_exp).^2))
    error_vector_2 = (mRNA_sigma28_exp .- mRNA_sigma28_sim)
    error_term_array[2,1] = 1.0*transpose(error_vector_2)*error_vector_2

    # mRNA C1 -
    #WM3 = diagm( 1.0 ./((mRNA_C1_std_exp).^2))
    error_vector_3 = (mRNA_C1_exp .- mRNA_C1_sim)
    error_term_array[3,1] = 1.0*transpose(error_vector_3)*error_vector_3

    # protein deGFP -
    #WM4 = diagm( 1.0 ./((mRNA_C1_std_exp).^2) )
    error_vector_4 = (protein_GFP_exp .- protein_GFP_sim)
    error_term_array[4,1] = 1.0*transpose(error_vector_4)*error_vector_4

    # the fifth term is an upper bound on C1 -
    prot_C1 = XSIM[:,9]
    UB = 100.0
    MV_C1 = maximum(prot_C1)
    error_term_C1 = 100000*max(0,(MV_C1 - UB))
    error_term_array[5,1] = error_term_C1

    # # 6th term - limit sigma28 mRNA peak to 200 -
    # mRNA_sigma28_actual = (1000)*XSIM[:,7]
    # UB_sigma28 = 800.0
    # MV_sigma28 = maximum(mRNA_sigma28_actual)
    # error_term_mRNA_sigma28_overshoot = 100000*max(0, (MV_sigma28 - UB_sigma28))
    # error_term_array[6,1] = error_term_mRNA_sigma28_overshoot
    # ===================================================================================================== #

    # error total -
    #error_total = sum(error_term_array)

    # return -
    return error_term_array
end

function main(path_to_data_dir::String, initial_parameter_array::Array{Float64,1}; rank_cutoff::Int64=4, maximum_number_of_iterations::Int64=100)

    # load the default data_dictionary -
    time_start = 0.0
    time_stop = 16.0
    time_step_size = 0.01

    # what is the host_type?
    host_type = :cell_free
    number_of_parameters = length(initial_parameter_array)

    # path to parameters -
    path_to_biophysical_constants_file = "./CellFree.json"

    # Load the data dictionary (uses the default biophysical_constants file)
    default_data_dictionary = build_data_dictionary((time_start,time_stop,time_step_size), path_to_biophysical_constants_file, host_type)

    # Set some values for the weights, gene lengths etc -
    model_data_dictionary = customize_data_dictionary(default_data_dictionary,host_type)

    # load the experimental data -
    exp_data_dictionary = load_experimental_data_dictionary(path_to_data_dir)

    # setup the functions -
    OF(P) = objective_function(P,time_start,time_step_size,time_stop,model_data_dictionary,exp_data_dictionary)
    NF(P) = neighbor_function(P;sigma=0.01)

    # make call to POETs -
    (EC,PC,RA) = estimate_ensemble(OF,NF,acceptance_probability_function,cooling_function,initial_parameter_array;rank_cutoff=rank_cutoff,maximum_number_of_iterations=maximum_number_of_iterations)

    # return -
    return (EC,PC,RA)
end

# setup -
path_to_data_dir = "$(pwd())/data"
pV = readdlm("pvec_initial_v5.dat")
pV = neighbor_function(pV; sigma=0.25)
EC = 0
PC = 0
RA = 0

# execute -
number_of_trials = 20
for trial_index = 1:number_of_trials

    global pV
    global EC
    global PC
    global RA

   # do a local step -
    if (mod(trial_index,2) == 0)

        # find the lowest score pV -
        sum_error_array = sum(EC,dims=1)
        best_p_index = argmin(vec(sum_error_array))
        pV_best = PC[:,best_p_index]

        # local refine -
        pV = local_refienment_step(path_to_data_dir, pV_best; iteration_max=500)
    end

    # main -
    (EC,PC,RA) = main(path_to_data_dir, vec(pV); rank_cutoff=3,maximum_number_of_iterations=100)

    # dump results to disk -
    fname = "./poets_ensemble_c1_restriction/RA_T$(trial_index).dat"
    writedlm(fname,RA)
    fname = "./poets_ensemble_c1_restriction/EC_T$(trial_index).dat"
    writedlm(fname,EC)
    fname = "./poets_ensemble_c1_restriction/PC_T$(trial_index).dat"
    writedlm(fname,PC)

    # find next (this will be overwritten by local search when index is even) -
    ECS = sum(EC,dims=1)
    idx_min = argmin(vec(ECS))
    pV = PC[:,idx_min]
end

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
"""
The barrel vault (Scordelis-Lo) roof is one of the benchmarks for linear elastic
analysis of shells.

The candidate element's usefulness in irregular geometries (and most practical
cases involve a high degree of geometric irregularity) is tested. As would be
expected,the irregular mesh results are not as good as those provided by a
regular meshwith the same number of variables.

Problem description

The physical basis of the problem is a deeply arched roof supported only
bydiaphragms at its curved edges (an aircraft hanger), deforming under its own
weight. It is interesting to observe that the geometry is such that the
centerpoint of the roof moves upward under the self-weight(downwardly directed)
load. Perhaps this is one reason why the problem is not straightforward
numerically.
"""
module scordelis_lo_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

# analytical solution for the vertical deflection and the midpoint of the
# free edge
analyt_sol = -0.3024
# Parameters:
E = 4.32e8
nu = 0.0
thickness = 0.25 # geometrical dimensions are in feet
R = 25.0
L = 50.0

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ)
    r[2] = 0.0
    r[3] += R
    csmatout[:, 3] .= vec(r) / norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function _execute_dsg_model(formul, n = 8, visualize = true)

    tolerance = R / n / 10
    # fens, fes = T3blockrand(40/360*2*pi,L/2,n,n);
    fens, fes = T3block(40 / 360 * 2 * pi, L / 2, n, n, :b)
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a = fens.xyz[i, 1]
        y = fens.xyz[i, 2]
        fens.xyz[i, :] .= (R * sin(a), y, R * (cos(a) - 1))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.drilling_stiffness_scale = 1.0e-4
    femm.mult_el_size = 0.2
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz, 1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz, 1), 6))

    # Apply EBC's
    # rigid diaphragm
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [1, 3, 5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf L / 2 L / 2 -Inf Inf], inflate = tolerance)
    for i in [2, 4, 6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1, 5, 6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi)

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi)

    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[sin(40 / 360 * 2 * pi) * 25 sin(40 / 360 * 2 * pi) * 25 L / 2 L / 2 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt[0, 0, -90, 0, 0, 0])
    F = distribloads(lfemm, geom0, dchi, fi, 3)

    # Solve
    U = K \ F
    scattersysvec!(dchi, U[:])
    result = dchi.values[nl, 3][1]
    @info "Solution: $(result), $(round(result/analyt_sol*100, digits = 4))%"

    # Visualization
    if visualize

        # Generate a graphical display of resultants
        scalars = []
        for nc in 1:3
            fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            push!(scalars, ("m$nc", fld.values))
            fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            push!(scalars, ("em$nc", fld.values))
        end
        vtkwrite("scordelis_lo_examples-$(n)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
        scalars = []
        for nc in 1:3
            fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            push!(scalars, ("n$nc", fld.values))
            fld = elemfieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            push!(scalars, ("en$nc", fld.values))
        end
        vtkwrite("scordelis_lo_examples-$(n)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
        scalars = []
        for nc in 1:2
            fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            push!(scalars, ("q$nc", fld.values))
            fld = elemfieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            push!(scalars, ("eq$nc", fld.values))
        end
        vtkwrite("scordelis_lo_examples-$(n)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

        vtkwrite("scordelis_lo_examples-$(n)-uur.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3]), ("ur", dchi.values[:, 4:6])])

        scattersysvec!(dchi, (L / 8) / maximum(abs.(U)) .* U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -L / 2]; [L / 2 L / 2 L / 2]]),
            #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end

    result
end

function test_convergence(ns = [4, 6, 8, 16, 32], visualize = false)
    formul = FEMMShellT3FFModule
    @info "Scordelis-Lo shell"
    results = []
    for n in ns
        v = _execute_dsg_model(formul, n, visualize)
        push!(results, v / (-0.3024) * 100)
    end
    return ns, results
end

end # module

using .scordelis_lo_examples

# Visualized internal resultants
ns, results = scordelis_lo_examples.test_convergence([2 * 16,], true)

# These results come from Table 9 of An efficient three‑node triangular
# Mindlin–Reissner flat shell element, Hosein Sangtarash1 · Hamed Ghohani
# Arab1 · Mohammad R. Sohrabi1 · Mohammad R. Ghasemi1
# they are all normalized relative to 0.3024# .   #
# Mesh Subdivision
# 4, 8, 16, 32

Allman = [
    1.004
    0.987
    0.987
    0.988] .* 100
Cook = [
    0.907
    0.929
    0.950
    0.981] .* 100
Providas_Kattis = [
    0.734
    0.815
    0.873
    0.967] .* 100
Shin_Lee = [
    1.379
    1.023
    1.004
    missing] .* 100
MITC3plus = [
    0.669
    missing
    0.857
    0.955] .* 100
TMRFS = [
    0.924
    0.963
    0.974
    0.998
] .* 100

using PGFPlotsX

objects = []

ns, results = scordelis_lo_examples.test_convergence()

compensate(r) = (0.3024 / 0.3006) .* r
all_results = [("Present", compensate(results), "*"), ("Allman", compensate(Allman), "x"), ("ProvKat", compensate(Providas_Kattis), "triangle"), ("Cook", compensate(Cook), "square"), ("SL", compensate(Shin_Lee), "o"), ("MITC3+", compensate(MITC3plus), "diamond")]

for r in all_results
    @pgf p = PGFPlotsX.Plot(
        {
            color = "black",
            line_width = 0.7,
            style = "solid",
            mark = "$(r[3])"
        },
        Coordinates([v for v in zip(ns, r[2]) if v[2] !== missing])
    )
    push!(objects, p)
    push!(objects, LegendEntry("$(r[1])"))
end

@pgf ax = Axis(
    {
        xlabel = "Number of Elements / side [ND]",
        ylabel = "Normalized Displacement [ND]",
        # xmin = range[1],
        # xmax = range[2],
        xmode = "linear",
        ymode = "linear",
        yminorgrids = "true",
        grid = "both",
        legend_pos = "north east"
    },
    objects...
)

display(ax)
pgfsave("scordelis_lo_examples-convergence.pdf", ax)

# ns, results = scordelis_lo_examples.test_convergence([4, 8, 16])
# q1, q2, q3 = results
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

# ns, results = scordelis_lo_examples.test_convergence(2 .* [4, 8, 16])
# q1, q2, q3 = results
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

ns, results = scordelis_lo_examples.test_convergence([16, 32, 64])
q1, q2, q3 = results
@show qtrue = (q2^2 - q1 * q3) / (2 * q2 - q1 - q3) * 0.3024 / 100

# ns, results = scordelis_lo_examples.test_convergence(8 .* [4, 8, 16])
# q1, q2, q3 = results
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

# ns, results = scordelis_lo_examples.test_convergence(16 .* [4, 8, 16])
# q1, q2, q3 = results
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

# ns, results = scordelis_lo_examples.test_convergence(32 .* [4, 8, 16])
# q1, q2, q3 = results
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

# ns, results = scordelis_lo_examples.test_convergence(64 .* [4, 8, 16])
# q1, q2, q3 = results
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

# julia> q1, q2, q3 = results
# 3-element Vector{Any}:
#  99.65925843568742
#  99.6820694134264
#  99.68801945695984

q1 = -0.2885517078616122
q2 = -0.2978284489578761
q3 = -0.3003728234489632
@show qtrue = (q2^2 - q1 * q3) / (2 * q2 - q1 - q3)

# julia> @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)
# qtrue = (q2 ^ 2 - q1 * q3) / ((2q2 - q1) - q3) = 99.69011916380857
# 99.69011916380857

# using FinEtools.AlgoBaseModule
# AlgoBaseModule.richextrapol(Float64.(results) , 1 ./ ns)
# (99.69011916375054, 1.9387561566847429, 1440.1118605939823, [3.750472155061857e-14, 3.239596091386687e-14, -1.0034074265918846e-14])

# # Data from Kiendl
# # Number of knots
# #    4                   8                   16                  32                  64
# 2  0.069489140101432   0.238098983889902   0.295702595727846   0.300237752650471   0.300558363037311
# 3  0.279944686988951   0.300065498003107   0.300584157108943   0.300592331224518   0.300592454396341
# 4  0.300000190520180   0.300589846259349   0.300592432120400   0.300592456541168   0.300592456632799
# 5  0.300966825455877   0.300589733200638   0.300592456543583   0.300592456635101   0.300592456617171

# # Suitable for extrapolation: p=5, first three results
# # Surprising: convergence from above?
# q1, q2, q3 = 0.300966825455877,   0.300589733200638,   0.300592456543583
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)  # = 0.30059243701673916
# # Not suitable for extrapolation: p=5, last three results
# q1, q2, q3 = 0.300592456543583,   0.300592456635101,   0.300592456617171
# @show qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)

# Abaqus, STRI3 elements, at regular mesh 64, 128, 256
q1 = 0.2997
q2 = 0.3003
q3 = 0.3005
@show qtrue = (q2^2 - q1 * q3) / (2 * q2 - q1 - q3) # 0.3006

# Abaqus, S4 elements, regular mesh 64, 128, 256
q1 = 0.3018
q2 = 0.301976
q3 = 0.302021
@show qtrue = (q2^2 - q1 * q3) / (2 * q2 - q1 - q3) # 0.302036

# T3FF, mesh with 128, 256, 512 elements per side
# julia> ns, results = scordelis_lo_examples.test_convergence([16, 32, 64].*8)
# [ Info: Scordelis-Lo shell
# [ Info: Solution: -3.01162e-01, 9.95907e+01%
# [ Info: Solution: -3.01378e-01, 9.96622e+01%
# [ Info: Solution: -3.01441e-01, 9.96828e+01%
q1, q2, q3 = [9.95907e+01, 9.96622e+01, 9.96828e+01]

qtrue = (((q2 ^ 2 - q1 * q3) / ((2q2 - q1) - q3)) * 3.02400e-01) / 100 # = 3.01466e-01


text:
using Test

using MIToS.Utils
using MIToS.MSA
using MIToS.Information
using MIToS.PDB
using MIToS.SIFTS
using MIToS.Pfam
using LinearAlgebra
using Random
using Statistics            # mean
using DelimitedFiles        # readdlm
using ROCAnalysis           # AUC
using Clustering            # test/MSA/Hobohm.jl
using NamedArrays           # array
using StatsBase             # WeightVec
using PairwiseListMatrices  # getlist

const DATA = joinpath(pwd(), "data")

# Utils
@testset "Utils" begin
    include(joinpath("Utils", "GeneralUtils.jl"))
end

# MSA
@testset "MSA" begin
    include(joinpath("MSA", "Residues.jl"))
    include(joinpath("MSA", "Alphabet.jl"))
    include(joinpath("MSA", "ThreeLetters.jl"))
    include(joinpath("MSA", "Annotations.jl"))
    include(joinpath("MSA", "MultipleSequenceAlignment.jl"))
    include(joinpath("MSA", "GeneralParserMethods.jl"))
    include(joinpath("MSA", "IO.jl"))
    include(joinpath("MSA", "General.jl"))
    include(joinpath("MSA", "MSAEditing.jl"))
    include(joinpath("MSA", "MSAStats.jl"))
    include(joinpath("MSA", "Shuffle.jl"))
    include(joinpath("MSA", "Identity.jl"))
    include(joinpath("MSA", "Hobohm.jl"))
    include(joinpath("MSA", "MSAAnnotations.jl"))
end

# Information
@testset "Information" begin
    include(joinpath("Information", "ContingencyTables.jl"))
    include(joinpath("Information", "Counters.jl"))
    include(joinpath("Information", "InformationMeasures.jl"))
    include(joinpath("Information", "Iterations.jl"))
    include(joinpath("Information", "CorrectedMutualInformation.jl"))
    include(joinpath("Information", "Gaps.jl"))
    include(joinpath("Information", "Externals.jl"))
end

# PDB
@testset "PDB" begin
    include(joinpath("PDB", "PDB.jl"))
    include(joinpath("PDB", "Contacts.jl"))
    include(joinpath("PDB", "Kabsch.jl"))
    include(joinpath("PDB", "Internals.jl"))
end

# SIFTS
@testset "SIFTS" begin
    include(joinpath("SIFTS", "SIFTS.jl"))
end

# Pfam
@testset "Pfam" begin
    include(joinpath("Pfam", "Pfam.jl"))
end

# Scripts
@testset "Scripts" begin
    include(joinpath("Scripts", "Template.jl"))
    include(joinpath("Scripts", "Scripts.jl"))
end

print("""

----- =D -----

""")


text:
struct KDEException <: Exception
    n::Int
    d::Int
    cov_det::Float64
    unique_values::Int
    KDEException(n, d, cov_det, unique_values) = new(n, d, cov_det, unique_values)
    KDEException(data) = new(size(data, 2), size(data, 1), det(cov(data; dims=2)), size(unique(data; dims=2), 2))
end

function Base.showerror(io::IO, kdee::KDEException)
    print(io, "KDE failed. ")
    if kdee.n < kdee.d
        print(io,"Cannot estimate covariance matrix if number of observations is less than the number of attributes." *
        "See https://en.wikipedia.org/wiki/Estimation_of_covariance_matrices. ")
    end
    if kdee.cov_det ≈ 0
        print(io,"""Covariance matrix is singular. """)
    end
    print(io, "Debug info: the number of attributes is $(kdee.d) and the number of unique observations is $(kdee.unique_values). ")
    return nothing
end

struct MissingLabelTypeException <: Exception
    t::Symbol
end

function Base.showerror(io::IO, e::MissingLabelTypeException)
    print(io, "QS failed. ")
    print(io, "Debug info: labels do not contain '$(e.t)'.")
    return nothing
end

function filter_array(x::Array{T, 2}, remove_indices::Array{Int, 1})::Array{T, 2} where T <: Real
    return x[:, setdiff(1:size(x, 2), remove_indices)]
end

function multi_kde(x::Array{Int, 2}, bw_method="scott")::PyCall.PyObject
    return multi_kde(convert.(Float64, x), bw_method)
end

function multi_kde(x::Array{Float64, 2}, bw_method="scott")::PyCall.PyObject
    size(x, 2) >= size(x, 1) || throw(KDEException(x))
    !(det(cov(x; dims=2)) ≈ 0) || throw(KDEException(x))
    return gaussian_kde(x, bw_method)
end

function leave_out_one_cv_kde(x::Array{T, 2}, bw_method="scott")::Float64 where T <: Real
    mean(multi_kde(x[:, 1:end .!= i], bw_method)(x[:, i]) for i in 1:size(x, 2))[1]
end

function knn_indices(x::Array{T, 2}; k=1::Int)::Array{Int, 2} where T <: Real
    kdtree = KDTree(x)
    return hcat(knn(kdtree, x, k + 1, true)[1]...)[2:end,:]
end

function knn_mean_dist(x::Array{T, 2}; k=1)::Array{Float64,1} where T <: Real
    kdtree = KDTree(x)
    return map(d -> mean(d[2:end]), knn(kdtree, x, k + 1, true)[2])
end

function initialize_qs(qs::DataType, model::OCClassifier, data::Array{T, 2}, params::Dict{Symbol, <:Any})::QueryStrategy where T <: Real
    if qs <: HybridQs
        return qs(model, data; params...)
    elseif qs <: ModelBasedQs
        return qs(model; params...)
    elseif qs <: DataBasedQs
        kernel = get_kernel(model)
        if typeof(kernel) == GaussianKernel
            return qs(data, bw_method=MLKernels.getvalue(strategy.kernel.alpha); params...)
        else
            return qs(data; params...)
        end
    elseif qs <: QueryStrategy
        return qs(; params...)
    end
    throw(ErrorException("Unknown query strategy of type $(qs)."))
end


text:
module TablesCatTest

using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#HPAT.DistributedPass.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.setCreateMain(true)
#ParallelAccelerator.CGen.set_debug_level(3)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)
#CompilerTools.CFGs.set_debug_level(3)

@acc hpat function test1(file_name)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    t2 = aggregate(t1,:uid = :userid, :ss=sum(:val2))
    t3 = aggregate(t1, :uid = :userid, :ss=sum(:val2>2.0))
    data = transpose(typed_hcat(Float64, t3[:uid], t3[:ss], t2[:ss]))
    return sum(data)
end

using HDF5
using MPI
using Base.Test

function main()

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        h5write("test1_1.hdf5", "/userid", [1, 2, 3, 1, 2])
        h5write("test1_1.hdf5", "/val2", [1.1, 2.1, 3.1, 3.2, 1.9])
    end
    s = test1("test1_1.hdf5")
    println(s)

    @test_approx_eq s 20.4

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        rm("test1_1.hdf5")
    end

end

end

TablesCatTest.main()


text:
#
# Either
#

"""
    Either{L,R}

A binary tagged union of `L` and `R`; either  [`Left`](@ref) or [`Right`](@ref).
"""
abstract type Either{L,R} end

for (E, c, V, O) in [(:Left, :left, :L, :R), (:Right, :right, :R, :L)]
    ise = Symbol(:is, c)

    # string and expr interpolations don't mix well
    dE = """
        $E{L,R} <: Either{L,R}

    $E case of `Either{L,R}`, constructed with [`$c`](@ref).
    """

    dc = """
        $c(v::$V, [$O])

    Return a `$E{L,R}(v) <: Either{L,R}`, with `$O = Any` as the default.
    """

    dis = """
        $ise(e)

    Return whether `e` is of type `$E`.
    """

    @eval begin
        @doc $dE
        struct $E{L,R} <: Either{L,R}
            v::$V
        end
        $E{$O}(v::$V) where {L,R} = $E{L,R}(v)

        Base.promote_rule(::Type{$E{L1,R1}}, ::Type{$E{L2,R2}}) where {L1,L2,R1,R2} =
                $E{promote_type(L1, L2), promote_type(R1, R2)}
        Base.convert(::Type{$E{L,R}}, e::$E) where {L,R} = $E{L,R}(value(e))

        @doc $dc
        $c(v, $O::Type=Any) = $E{$O}(v)
        $c(v, ::$O) where {$O} = $E{$O}(v)

        @doc $dis
        $ise(::$E) = true
        $ise(_) = false
    end
end

lefttype(::Either{L}) where {L} = L
righttype(::Either{L,R}) where {L,R} = R

Base.eltype(l::Left) = lefttype(l)
Base.eltype(r::Right) = righttype(r)

value(e::Either) = e.v

Base.join(e::Either) = _join(e, value(e))
_join(e, v::Either) = _join(v, value(v))
_join(e, v) = e

"""
    map(f, e::Either, [R2])

Return `right(f(value(e)))` if `e::Right`, else return `e`

`R2` is `f`s expected return type, to ensure the map returns `e::Left{L, R2}`.
"""
Base.map(f, e::Either) = Base.map(f, e, righttype(e))
Base.map(_, l::Left, R::Type) = left(l, R)
Base.map(f, r::Right, _::Type) = right(f(value(r)), lefttype(r))

Base.show(io::IO, l::Left) =
        print(io, "Left($(value(l)), ::$(righttype(l)))")
Base.show(io::IO, r::Right) =
        print(io, "Right(::$(lefttype(r)), $(value(r)))")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

const name = "libcgal-julia"

version = v"0.99.0"
const matches = match(r"refs/tags/(.*)", get(ENV, "GITHUB_REF", ""))
if matches !== nothing
    try
        global version = VersionNumber(first(matches.captures))
    catch e
        @warn e.msg
    end
end
@info "$name v$version"

# Collection of sources required to build CGAL
const sources = [
     "$(@__DIR__)/../../"
]

# Dependencies that must be installed before this package can be built
const dependencies = [
    "CGAL_jll",
    "libcxxwrap_julia_jll",
]

# Bash recipe for building across all platforms
const script = raw"""
## pre-build setup
# exit on error
set -eu

# HACK: download julia..
curl -Lo julia.tar.gz https://github.com/JuliaPackaging/JuliaBuilder/releases/download/v1.0.0-2/julia-1.0.0-$target.tar.gz
mkdir julia && tar xf julia.tar.gz -C julia
Julia_PREFIX="$PWD/julia"

## configure build
mkdir -p build && cd build

cmake ../ \
  `# cmake specific` \
  -DCMAKE_TOOLCHAIN_FILE="$CMAKE_TARGET_TOOLCHAIN" \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_CXX_FLAGS="-march=x86-64" \
  -DCMAKE_FIND_ROOT_PATH="$prefix" \
  -DCMAKE_INSTALL_PREFIX="$prefix" \
  `# tell libcxxwrap-julia where julia is` \
  -DJulia_PREFIX="$Julia_PREFIX"

## and away we go..
VERBOSE=ON cmake --build . --config Release --target install -- -j$nproc

install_license ../LICENSE

# HACK: Apparently, this isn't a simple build system anymore..
case $target in
  *mingw32*) mv "$prefix/lib/"*.dll "$prefix/bin" ;;
esac
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
const platforms = [
    Linux(:x86_64, libc=:glibc),
    Windows(:i686),
    Windows(:x86_64),
    MacOS(:x86_64),
] |> expand_cxxstring_abis

# The products that we will ensure are always built
const products = [
    LibraryProduct("libcgal_julia_exact", :libcgal_julia_exact),
    LibraryProduct("libcgal_julia_inexact", :libcgal_julia_inexact),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: rgcv/libcgal-julia

max_issues_repo_name: rgcv/libcgal-julia

max_forks_repo_name: rgcv/libcgal-julia


text:
"""
    mese!(spin, TR, TE, nechoes; αex, αinv)

Simulate the steady-state signal acquired from a multi-echo spin echo (MESE)
sequence, assuming instantaneous excitations and ideal spoiling.

# Arguments
- `spin::AbstractSpin`: Spin whose signal to acquire
- `TR::Real`: Repetition time of the sequence (ms)
- `TE::Real`: Echo time of the sequence (ms)
- `nechoes::Integer`: Number of spin echoes to acquire
- `αex::Real = π/2`: Flip angle for the initial excitation (rad)
- `αinv::Real = π`: Flip angle for the refocussing (inversion) pulse (rad)

# Return
- `signal::Vector{ComplexF64}`: Steady-state signal at each of the `nechoes`
    echo times
"""
function mese!(spin::AbstractSpin, TR::Real, TE::Real, nechoes::Integer;
               αex::Real = π/2, αinv::Real = π)

    TR >= TE * nechoes ||
        throw(ArgumentError("TR must be greater than or equal to TE * nechoes"))

    # Precompute spin dynamics
    Dex = excitation(spin, 0, αex)
    Dinv = excitation(spin, -π/2, αinv)
    Dte = freeprecess(spin, TE/2)
    Dtr = freeprecess(spin, TR - TE * nechoes)
    Decho = combine(Dte, Dinv, Dte)
    Dspoil = (spoil(spin), zeros(length(spin.M)))

    # Calculate steady-state magnetization immediately following excitation
    D = Decho
    for e = 2:nechoes
        D = combine(D, Decho)
    end
    (A, B) = combine(D, Dtr, Dspoil, Dex)
    spin.M[:] = (Diagonal(ones(Bool, size(A, 1))) - A) \ B

    # Calculate steady-state signal at each echo
    signal = zeros(ComplexF64, nechoes)
    for e = 1:nechoes
        applydynamics!(spin, Decho...)
        signal[e] = spin.signal
    end

    return signal

end

"""
    mese(spin, TR, TE, nechoes; αex, αinv, nspins, ncycles)

Simulate the steady-state signal acquired from a multi-echo spin echo (MESE)
sequence, assuming instantaneous excitations and nonideal gradient spoiling.

The output is a vector of `nechoes` signal values that are the complex mean of
`nspins` copies of `spin` at slightly offset positions to induce `ncycles`
cycles of phase across the spins.

# Arguments
- `spin::AbstractSpin`: Spin whose signal to acquire
- `TR::Real`: Repetition time of the sequence (ms)
- `TE::Real`: Echo time of the sequence (ms)
- `nechoes::Integer`: Number of spin echoes to acquire
- `αex::Real = π/2`: Flip angle for the initial excitation (rad)
- `αinv::Real = π`: Flip angle for the refocussing (inversion) pulse (rad)
- `nspins::Integer = 10`: Number of intra-voxel spins to simulate
- `ncycles::Real = 1`: Number of cycles of phase across a voxel induced by each
    spoiler gradient (typically integer-valued)

# Return
- `signal::Vector{ComplexF64}`: Ensemble steady-state signal at each of the
    `nechoes` echo times
"""
function mese(spin::AbstractSpin, TR::Real, TE::Real, nechoes::Integer;
              αex::Real = π/2, αinv::Real = π, nspins::Integer = 10,
              ncycles::Real = 1)

    # Pick an arbitrary gradient strength and duration and compute the spatial
    # locations that will provide a uniform ncycles of phase
    gradz = 0.3 # G/cm
    Tg = TE / 4 # ms
    zmax = ncycles * 2π / (GAMMA * gradz * Tg/1000) # cm
    z = (1:nspins)/nspins  * zmax

    # Make nspins copies of the provided spin at different spatial locations
    spins = varyz(spin, z)

    # Compute the MESE signal for each spin
    signals = map(spin -> mese!(spin, TR, TE, nechoes, [0,0,gradz], Tg,
                            αex = αex, αinv = αinv), spins) # [nspins][nechoes]

    # Find the average signal for each echo
    signal = sum(signals) / nspins

    return signal

end

"""
    varyz(spin, z)

Copy `spin`, offsetting the z position of each copy by elements of z.
"""
function varyz(spin::Spin, z)

    map(z -> Spin(spin.M0, spin.T1, spin.T2, spin.Δf, spin.pos + [0,0,z]), z)

end

function varyz(spin::SpinMC, z)

    map(z -> SpinMC(spin.M0, spin.frac, spin.T1, spin.T2, spin.Δf, spin.τ,
                    spin.pos + [0,0,z]), z)

end

"""
    mese!(spin, TR, TE, nechoes, grad, Tg; αex, αinv)

Simulate the steady-state signal acquired from a multi-echo spin echo (MESE)
sequence, assuming instantaneous excitations and nonideal gradient spoiling.

The spoiler gradient is played before and after each refocusing pulse. There is
no RF spoiling. The phase of the excitation pulse is not reversed after each TR.

# Arguments
- `spin::AbstractSpin`: Spin whose signal to acquire
- `TR::Real`: Repetition time of the sequence (ms)
- `TE::Real`: Echo time of the sequence (ms)
- `nechoes::Integer`: Number of spin echoes to acquire
- `grad::AbstractArray{<:Real}`: Spoiler gradient amplitudes [gx, gy, gz] (G/cm)
- `Tg::Real`: Spoiler gradient duration (ms)
- `αex::Real = π/2`: Flip angle for the initial excitation (rad)
- `αinv::Real = π`: Flip angle for the refocussing (inversion) pulse (rad)

# Return
- `signal::Vector{ComplexF64}`: Steady-state signal at each of the `nechoes`
    echo times
"""
function mese!(spin::AbstractSpin, TR::Real, TE::Real, nechoes::Integer,
               grad::AbstractArray{<:Real}, Tg::Real; αex::Real = π/2,
               αinv::Real = π)

    TR >= TE * nechoes ||
        throw(ArgumentError("TR must be greater than or equal to TE * nechoes"))
    Tg <= TE / 2 ||
        throw(ArgumentError("Tg must be less than or equal to TE / 2"))

    # Precompute spin dynamics
    Dex = excitation(spin, 0, αex)
    Dinv = excitation(spin, -π/2, αinv)
    Dg = freeprecess(spin, Tg, grad)
    Dte = freeprecess(spin, TE/2 - Tg)
    Dtr = freeprecess(spin, TR - TE * nechoes)
    Decho = combine(Dte, Dg, Dinv, Dg, Dte)

    # Calculate steady-state magnetization immediately following excitation
    D = Decho
    for e = 2:nechoes
        D = combine(D, Decho)
    end
    (A, B) = combine(D, Dtr, Dex)
    spin.M[:] = (Diagonal(ones(Bool, size(A, 1))) - A) \ B

    # Calculate steady-state signal at each echo
    signal = zeros(ComplexF64, nechoes)
    for e = 1:nechoes
        applydynamics!(spin, Decho...)
        signal[e] = spin.signal
    end

    return signal

end

"""
    spgr!(spin, TR, TE, α)

Simulate the steady-state signal acquired from a spoiled gradient echo (SPGR)
sequence, assuming instantaneous excitations and ideal spoiling.

# Arguments
- `spin::AbstractSpin`: Spin whose signal to acquire
- `TR::Real`: Repetition time of the sequence (ms)
- `TE::Real`: Echo time of the sequence (ms)
- `α::Real`: Flip angle (rad)

# Return
- `signal::ComplexF64`: Steady-state signal
"""
function spgr!(spin::AbstractSpin, TR::Real, TE::Real, α::Real)

    TR >= TE || throw(ArgumentError("TR must be greater than or equal to TE"))

    # Precompute spin dynamics
    Dex = excitation(spin, 0, α)
    Dte = freeprecess(spin, TE)
    Dtr = freeprecess(spin, TR - TE)
    Dspoil = (spoil(spin), zeros(length(spin.M)))

    # Calculate steady-state magnetization immediately following excitation
    (A, B) = combine(Dte, Dtr, Dspoil, Dex)
    spin.M[:] = (Diagonal(ones(Bool, size(A, 1))) - A) \ B

    # Calculate steady-state signal at echo time
    applydynamics!(spin, Dte...)

    return spin.signal

end

"""
    spgr!(spin, TR, TE, α, grad, Tg; Δθinc, nTR)

Simulate the steady-state signal acquired from a spoiled gradient echo (SPGR)
sequence, assuming instantaneous excitations and nonideal spoiling.

The steady-state signal is a true steady-state when `Δθinc` is 0, otherwise
`nTR` TR's are simulated to produce a pseudo steady-state. (Note, it is only a
pseudo steady-state when many spins are averaged together.)

# Arguments
- `spin::AbstractSpin`: Spin whose signal to acquire
- `TR::Real`: Repetition time of the sequence (ms)
- `TE::Real`: Echo time of the sequence (ms)
- `α::Real`: Flip angle (rad)
- `grad::AbstractArray{<:Real}`: Spoiler gradient amplitudes [gx, gy, gz] (G/cm)
- `Tg::Real`: Spoiler gradient duration (ms)
- `Δθinc::Real = deg2rad(117)`: RF phase increment for RF spoiling (rad)
- `nTR::Integer = 100`: Number of TR's to simulate when `Δθinc` is not 0

# Return
- `signal::ComplexF64`: Steady-state signal
"""
function spgr!(spin::AbstractSpin, TR::Real, TE::Real, α::Real,
               grad::AbstractArray{<:Real}, Tg::Real;
               Δθinc::Real = deg2rad(117), nTR::Integer = 100)

    # I would prefer to split this into two functions, but it is not possible
    # (without having multiple names like spgr_rfspoil) because I want Δθinc and
    # nTR to be optional arguments, but doing so would case the RF-spoiled case
    # and the non-RF-spoiled case to have the same method signature. And I like
    # the simplicity of having just one spgr function, so I don't want to add,
    # e.g., spgr_rfspoil.

    TR >= TE + Tg ||
        throw(ArgumentError("TR must be greater than or equal to TE + Tg"))

    # Precompute spin dynamics
    Dte = freeprecess(spin, TE)
    Dtr = freeprecess(spin, TR - TE - Tg)
    Dtg = freeprecess(spin, Tg, grad)
    D = combine(Dte, Dtr, Dtg)

    if Δθinc == 0

        # Precompute spin dynamics
        Dex = excitation(spin, 0, α)

        # Calculate steady-state magnetization immediately following excitation
        (A, B) = combine(D, Dex)
        spin.M[:] = (Diagonal(ones(Bool, size(A, 1))) - A) \ B

        # Calculate steady-state signal at echo time
        applydynamics!(spin, Dte...)

        return spin.signal

    else

        # Initialize RF spoiling parameters
        θ = 0
        Δθ = Δθinc

        # Simulate nTR TR's
        for rep = 1:nTR

            excitation!(spin, θ, α)
            applydynamics!(spin, D...)
            θ += Δθ
            Δθ += Δθinc

        end

        # Calculate signal at echo time
        excitation!(spin, θ, α)
        applydynamics!(spin, Dte...)

        return spin.signal * exp(im * θ)

    end

end


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 9d2617c5-ebf2-46c8-910f-190ba1d160ac
using BenchmarkTools

# ╔═╡ 89196053-24db-4013-a1e3-3d3b22c4eedc
using Combinatorics

# ╔═╡ 27547863-22e1-466e-951e-2a5677bb9501
using IterTools

# ╔═╡ bade5453-0522-4cbc-8200-78ec160f1192
using ConstraintSolver, JuMP

# ╔═╡ 568f057a-126d-11ec-23cb-f98bea4a534f
md"# Project Euler 1

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

"

# ╔═╡ daad198c-10c8-4e4c-93c4-d330ab22f3dd
naive1(x)=sum([i for i in 1:(x-1) if i%3==0 || i%5==0])

# ╔═╡ 2cca4797-d1ee-4c27-8243-736e3e468dd1
naive1(1000)

# ╔═╡ 6b82e993-2f3d-4d5c-8709-57050c7211de
sumRange(i,j) = (j-i+1)*(i+j)/2

# ╔═╡ 382a2b3a-4059-4cbe-8485-b0ecce99c451
sumRange(1,10)

# ╔═╡ e59e2718-acc6-4a98-8ec9-61cf0a375f90
md"We can also solve it with the inclusion-exclusion principle"

# ╔═╡ f7acd079-f79b-4967-a9ee-fc4965e2a06c
inclusionExclusionPrincipleProb1(n) = begin
	k=n-1
	sumRange(1,floor(k/3))*3+sumRange(1,floor(k/5))*5-sumRange(1,floor(k/15))*15
end

# ╔═╡ c118a05d-d1bb-4767-9f3a-0b5cefd69d47
inclusionExclusionPrincipleProb1(1000)

# ╔═╡ 078cb436-e238-4a69-bde6-156624db601b
@benchmark sum(1:1000000)

# ╔═╡ cea0b9f0-b908-4c33-8a5d-6a3019cfa47a
@benchmark sumRange(1,1000000)

# ╔═╡ 07d10397-f053-45e0-ab2b-6565c18db5a2
@benchmark inclusionExclusionPrincipleProb1(10000)

# ╔═╡ 464e87db-a12f-4e8b-a03a-ae1f943a42d0
@benchmark naive1(10000)

# ╔═╡ e323b98a-47a7-4182-9892-4d0b1891e295
md" so we can see we have about $(6.380*10^(-6)/(0.021*10^(-9))) speed-up"

# ╔═╡ 82baf161-7a98-4bb5-a928-b2273e4360cb
md"We can also do it for a general problem with the inclusion- exclusion principle as follows:
"

# ╔═╡ c7667675-8018-40f5-a1d9-8756f29466e0
prod([1,2,3])

# ╔═╡ ed4ded6b-925a-46e7-bcab-0cf0dc0fe80a
inclusionExclusionSumUpTo(n,multiples) = begin
	k=n-1
	res = 0
	for comb in combinations(multiples)
		pcomb=prod(comb)
		res += (-1)^(length(comb)-1)* sumRange(1,floor(k/pcomb))*pcomb
	end
	res
end

# ╔═╡ fa4d516c-74ea-44cb-bc18-e2b6cde80515
inclusionExclusionSumUpTo(1000,[3,5])

# ╔═╡ 120de079-e8e5-40b3-9582-c3741cf024cb
@benchmark inclusionExclusionSumUpTo(100000,[3,5,7,11,13])

# ╔═╡ 6ea414f8-d5bb-4432-9afb-6d2b2cfd1d01
md"We can see it's pretty fast..."

# ╔═╡ 6d793a48-b48e-4547-ac87-c23a1b03de9c
md" # Problem 2
Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

"

# ╔═╡ e0470dad-e1ab-4603-85d0-e165f80c2463
md"Copy some fibonacci channel code from [this blog post about lazy sequences](https://www.oxinabox.net/2017/11/18/Lazy-Sequences-in-Julia.html)"

# ╔═╡ c6a079d7-2773-41ba-81c8-773efaabd010
function fib_ch(buff=1)
    prev=0
    cur=1
    Channel(ctype=Int, csize=buff) do c
        while(true)
            put!(c, cur)
            prev, cur = (cur, prev+cur) # nice little avoiding a temp here
        end
    end
end

# ╔═╡ 61d320a0-2f4f-4cf4-b800-20bb77560c17

# ╔═╡ 2eb91ffc-799d-496e-8b50-6be2d497beaa
sum([x for x in collect(IterTools.takewhile(x->x<=4_000_000,fib_ch(100))) if x%2==0])

# ╔═╡ 5e201d20-3600-4b21-a764-8dd1e89cdf66
md" It's also quite fast!"

# ╔═╡ 6678078e-f30c-4e5e-8de5-f697d6c241fd
eulerProblem2WithLazySeq(limit)=sum([x for x in collect(IterTools.takewhile(x->x<=limit,fib_ch(100))) if x%2==0])

# ╔═╡ c2c72029-ffd2-463c-8548-bbe306178bd1
eulerProblem2WithLazySeq(4_000_000)

# ╔═╡ 77f86b91-0ac4-4ba4-a687-0ea49a0ddd10
@benchmark eulerProblem2WithLazySeq(4_000_000)

# ╔═╡ 314415de-ac90-479a-8665-fa9e51323740
@benchmark eulerProblem2WithLazySeq(1_000_000_000_000)

# ╔═╡ cea3d8d2-70e4-49c7-9698-d4085976e924
md"Let's figure out also how to solve some recurrent equations and do it like that

$a_n=a_{n-1}+a_{n-2}$

The characteristic equation of the generating function is

$x^2-x-1==0$, with roots

$x_1 = (1+√5)/2, x_2=(1-√5)/2$

Then we can search for a solution of the form $dx_1^n+ex_2^n$,
since then
$a_n-a_{n-1}-a_{n-2} = (dx_1^n+ex_2^n)-(dx_1^{n-1}+ex_2^{n-1})-(...) = dx_1^{n-2}(...=0)-ex_2^{n-2}(...=0)$
and setting $a_0 = 0, a_1=1$ and solving for
$d,e$ we get eventually Binet's formula:

$ϕ = (\sqrt(5)+1)/2$
$a_n = (ϕ^n-((-ϕ)^{-n}))/√5$

"

# ╔═╡ fba39823-d264-4210-80c8-491eb19486d2
md"# Product 3

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

"

# ╔═╡ d6a95c3b-30c1-4004-884c-2acdf4a2ef7f
# let's implement erathostens' sieve
function get_primes(n::Int64)
	numbers = fill(true,n)
	numbers[1] = false

	for i = 2:isqrt(n)
		if numbers[i]
			for j = i:div(n,i)
				numbers[i*j] = false
			end
		end
	end
	primes = findall(numbers)
	return primes
end

# ╔═╡ 080bd517-473a-4c13-9969-a4775a2acc16

# ╔═╡ 4746f780-750d-4151-bdfb-9687456a36f5
function findFactors(n::Int64)
 all_primes = get_primes(isqrt(600851475143))
 res = []
 for prime in all_primes
		if mod(n,prime)==0
			push!(res,prime)
			while mod(n,prime)==0
				n= (n//prime)
			end
		end
	end
	res
end

# ╔═╡ f2279703-313c-4146-8e01-52c5fecf8cfd
findFactors(600851475143)

# ╔═╡ 69853164-2ca0-4cba-9aef-307f0766b273
md"
# Problem 4

 The largest palindrome made fromt he product of 2 2-digit numbers is 9009=91*99.
What's the largest palindorme made forom the product of two 3-digit numbers...

"

# ╔═╡ 3c93ec46-633e-4d97-aa53-6d16b7e85502
CS = ConstraintSolver

# ╔═╡ 1c57967f-f38f-41e2-bc7c-51d1dd681e0b
function CSPalindrome()
	mod = Model(CS.Optimizer)
	@variable(mod,0<=a[i=1:3]<=9,Int)
	@variable(mod,0<=b[i=1:3]<=9,Int)
	@variable(mod,0<=res[i=1:6]<=9,Int)
	@variable(mod,0<=carry[i=1:5]<=9,Int)
	@constraint(mod,three_digits_1,a[1]>=1)
	@constraint(mod,three_digits_2,b[1]>=1)
	@constraint(mod,first_digit,res[1]==b[1]*a[1])
	@objective(mod,Max,1_000_00*res[1]+10_000*res[2]+1000*res[3]+100*res[4]+10*res[5]+res[6])
	mod

end

# ╔═╡ b5da5b91-de7c-4b69-98f5-737a8a9ced34
begin
	md = CSPalindrome()
	optimize!(md)
end

# ╔═╡ e0d8e8d1-8e92-4ea3-9b38-bfb7378fe004
objective_value(md)

# ╔═╡ f6b8b94d-0176-4b19-9c81-2959cfff60e7
isPalindrome(k,l)= begin
 str1= string(k*l)
  str1==reverse(str1)
end

# ╔═╡ cc6dbc26-db1f-4534-95b9-7a801139e1d0
findPalindrome()=

	for i = 999:-1:100
	for j in i:-1:100
		if isPalindrome(i,j)
			return i,j
	end

	end
	end

# ╔═╡ c5d9bd7f-1a75-482f-a2aa-f653fdce083c
findPalindrome()

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
ConstraintSolver = "e0e52ebd-5523-408d-9ca3-7641f1cd1405"
IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"

[compat]
BenchmarkTools = "~1.1.4"
Combinatorics = "~1.0.2"
ConstraintSolver = "~0.7.0"
IterTools = "~1.3.0"
JuMP = "~0.21.10"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-beta3.0"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

[[deps.Arpack]]
deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
version = "0.4.0"

[[deps.Arpack_jll]]
deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
uuid = "68821587-b530-5797-8361-c406ea357684"
version = "3.5.0+3"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
git-tree-sha1 = "42ac5e523869a84eac9669eaceed9e4aa0e1587b"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.1.4"

[[deps.Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "30ee06de5ff870b45c78f529a6b093b3323256a3"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.1"

[[deps.CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

[[deps.CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "6071cb87be6a444ac75fdbf51b8e7273808ce62f"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.35.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstraintProgrammingExtensions]]
deps = ["LinearAlgebra", "MathOptInterface", "Test"]
git-tree-sha1 = "989acd05e57112a5fd1bc8849ebc398a614aa010"
uuid = "b65d079e-ed98-51d9-b0db-edee61a5c5f8"
version = "0.3.2"

[[deps.ConstraintSolver]]
deps = ["ConstraintProgrammingExtensions", "DataStructures", "Formatting", "JSON", "JuMP", "LightGraphs", "MathOptInterface", "MatrixNetworks", "Random", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "47a4c07efa6890e7da286495e22742eb19f3eed9"
uuid = "e0e52ebd-5523-408d-9ca3-7641f1cd1405"
version = "0.7.0"

[[deps.DataAPI]]
git-tree-sha1 = "bec2532f8adb82005476c141ec23e921fc20971b"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.8.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[deps.DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "b5e930ac60b613ef3406da6d4f42c35d8dc51419"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.19"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "60ed5f1643927479f845b0135bb369b031b541fa"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.14"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[deps.IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.JSONSchema]]
deps = ["HTTP", "JSON", "URIs"]
git-tree-sha1 = "2f49f7f86762a0fbbeef84912265a1ae61c4ef80"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.4"

[[deps.JuMP]]
deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Printf", "Random", "SparseArrays", "SpecialFunctions", "Statistics"]
git-tree-sha1 = "4358b7cbf2db36596bdbbe3becc6b9d87e4eb8f5"
uuid = "4076af6c-e467-56ae-b986-b466b2749572"
version = "0.21.10"

[[deps.KahanSummation]]
deps = ["Test"]
git-tree-sha1 = "1f01068b28d3ad83d4d1212a0ce8d7ecacb33482"
uuid = "8e2b3108-d4c1-50be-a7a2-16352aec75c3"
version = "0.1.0"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "86197a8ecb06e222d66797b0c2d2f0cc7b69e42b"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.2"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "575644e3c05b258250bb599e57cf73bbf1062901"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.9.22"

[[deps.MatrixNetworks]]
deps = ["Arpack", "DataStructures", "DelimitedFiles", "IterTools", "KahanSummation", "LinearAlgebra", "Printf", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "f6794d987d426bab570aabb2d03e2f40ca4fc438"
uuid = "4f449596-a032-5618-b826-5a251cb6dc11"
version = "1.0.2"

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "2ca267b08821e86c5ef4376cffed98a46c2cb205"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[deps.NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "46d7ccc7104860c38b11966dd1f72ff042f382e4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.10"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═568f057a-126d-11ec-23cb-f98bea4a534f
# ╠═daad198c-10c8-4e4c-93c4-d330ab22f3dd
# ╠═2cca4797-d1ee-4c27-8243-736e3e468dd1
# ╠═6b82e993-2f3d-4d5c-8709-57050c7211de
# ╠═382a2b3a-4059-4cbe-8485-b0ecce99c451
# ╠═e59e2718-acc6-4a98-8ec9-61cf0a375f90
# ╠═f7acd079-f79b-4967-a9ee-fc4965e2a06c
# ╠═c118a05d-d1bb-4767-9f3a-0b5cefd69d47
# ╠═9d2617c5-ebf2-46c8-910f-190ba1d160ac
# ╠═078cb436-e238-4a69-bde6-156624db601b
# ╠═cea0b9f0-b908-4c33-8a5d-6a3019cfa47a
# ╠═07d10397-f053-45e0-ab2b-6565c18db5a2
# ╠═464e87db-a12f-4e8b-a03a-ae1f943a42d0
# ╠═e323b98a-47a7-4182-9892-4d0b1891e295
# ╠═82baf161-7a98-4bb5-a928-b2273e4360cb
# ╠═89196053-24db-4013-a1e3-3d3b22c4eedc
# ╠═c7667675-8018-40f5-a1d9-8756f29466e0
# ╠═ed4ded6b-925a-46e7-bcab-0cf0dc0fe80a
# ╠═fa4d516c-74ea-44cb-bc18-e2b6cde80515
# ╠═120de079-e8e5-40b3-9582-c3741cf024cb
# ╠═6ea414f8-d5bb-4432-9afb-6d2b2cfd1d01
# ╠═6d793a48-b48e-4547-ac87-c23a1b03de9c
# ╠═e0470dad-e1ab-4603-85d0-e165f80c2463
# ╠═c6a079d7-2773-41ba-81c8-773efaabd010
# ╠═61d320a0-2f4f-4cf4-b800-20bb77560c17
# ╠═27547863-22e1-466e-951e-2a5677bb9501
# ╠═2eb91ffc-799d-496e-8b50-6be2d497beaa
# ╠═5e201d20-3600-4b21-a764-8dd1e89cdf66
# ╠═6678078e-f30c-4e5e-8de5-f697d6c241fd
# ╠═c2c72029-ffd2-463c-8548-bbe306178bd1
# ╠═77f86b91-0ac4-4ba4-a687-0ea49a0ddd10
# ╠═314415de-ac90-479a-8665-fa9e51323740
# ╠═cea3d8d2-70e4-49c7-9698-d4085976e924
# ╠═fba39823-d264-4210-80c8-491eb19486d2
# ╠═d6a95c3b-30c1-4004-884c-2acdf4a2ef7f
# ╠═080bd517-473a-4c13-9969-a4775a2acc16
# ╠═4746f780-750d-4151-bdfb-9687456a36f5
# ╠═f2279703-313c-4146-8e01-52c5fecf8cfd
# ╠═69853164-2ca0-4cba-9aef-307f0766b273
# ╠═bade5453-0522-4cbc-8200-78ec160f1192
# ╠═3c93ec46-633e-4d97-aa53-6d16b7e85502
# ╠═1c57967f-f38f-41e2-bc7c-51d1dd681e0b
# ╠═b5da5b91-de7c-4b69-98f5-737a8a9ced34
# ╠═e0d8e8d1-8e92-4ea3-9b38-bfb7378fe004
# ╠═f6b8b94d-0176-4b19-9c81-2959cfff60e7
# ╠═cc6dbc26-db1f-4534-95b9-7a801139e1d0
# ╠═c5d9bd7f-1a75-482f-a2aa-f653fdce083c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
function trans!(state, r, a, b)
    x = state[1]
    y = state[2]
    x_out = x * exp( r - a*x - b*y )
    y_out = y * exp( r - a*y + b*x )

    [x_out, y_out]
end

function simulate(inits, r, a, b, nburn=Int(1e3), nkeep=Int(50e3))

    state = deepcopy(inits)

    for i = 1:nburn
        state = trans!(state, r, a, b)
    end
    println(state)

    sim = zeros(Float64, nkeep, 2)

    for i = 1:nkeep
        state = trans!(state, r, a, b)
        sim[i,:] = copy(state)
    end

    sim
end

using Random
Random.seed!(1)

## dynamics are deterministic
sim = simulate([5.0, 5.0], 2.75, 0.5, 0.07, 1000, 10010)

using Plotly

nplot = 1000

plot(sim[1:nplot,1])
plot(sim[1:nplot,1], sim[2:(nplot+1),1], mode="markers")

trace = scatter3d(Dict(
    [ :x => sim[1:(nplot),1], :y => sim[1:(nplot),2], :z => sim[2:(nplot+1),1],
      :mode => "markers",
      :marker => Dict(
        :size => 2
      ) ]
))

Plotly.plot( trace )

trace2 = scatter3d(Dict(
    [ :x => sim[1:(nplot),1], :y => sim[2:(nplot+1),1], :z => sim[3:(nplot+2),1],
      :mode => "markers",
      :marker => Dict(
        :size => 2
      ) ]
))

Plotly.plot( trace2 )

using BayesInference
x = embed(sim[:,1], 10)

y = deepcopy(x[:,1])
X = deepcopy(x[:,2:11])

using RCall
@rput y X
R"X = log(X); y = log(y); save(file='Ricker2d_log.rda', X, y)"

using BSON
bson("Ricker2d.bson", y=deepcopy(y), X=deepcopy(X))
bson("Ricker2d_log.bson", y=deepcopy(log.(y)), X=deepcopy(log.(X)))

### plots
using Plotly

r = 2.75
a = 0.5
b = 0.07

nplot = 3000

trace_sysdat = scatter3d(Dict(
    :x => (sim[1:nplot,1]),
    :y => (sim[1:nplot,2]),
    :z => (sim[2:(nplot+1),1]),
    :marker => Dict(
        :color => "rgba(0,0,0,0.70)",
        :size => 2
    ),
    :mode => "markers"
))

trace_logsysdat = scatter3d(Dict(
    :x => log.(sim[1:nplot,1]),
    :y => log.(sim[1:nplot,2]),
    :z => log.(sim[2:(nplot+1),1]),
    :marker => Dict(
        :color => "rgba(0,0,0,0.70)",
        :size => 2
    ),
    :mode => "markers"
))

trace_tdedat = scatter3d(Dict(
    :x => (X[1:nplot,1]),
    :y => (X[1:nplot,2]),
    :z => (y[1:nplot]),
    :marker => Dict(
        :color => "rgba(0,0,0,0.70)",
        :size => 2
    ),
    :mode => "markers"
))

trace_logtdedat = scatter3d(Dict(
    :x => log.(X[1:nplot,1]),
    :y => log.(X[1:nplot,2]),
    :z => log.(y[1:nplot]),
    :marker => Dict(
        :color => "rgba(0,0,0,0.70)",
        :size => 2
    ),
    :mode => "markers"
))

function sys_surf(x, y, r, a, b)
    x_out = x * exp( r - a*x - b*y )
    y_out = y * exp( r - a*y + b*x )

    [x_out, y_out]
end

function sys_surflog(lx, ly, r, a, b)
    x = exp(lx)
    y = exp(ly)
    lx_out = lx + r - a*x - b*y
    ly_out = ly + r - a*y + b*x

    [lx_out, ly_out]
end

function tde_surf(x1, x2, r, a, b)
    lx1 = log(x1)
    lx2 = log(x2)
    ax2rl = a*x2 - r + lx1 - lx2
    x1 * exp( r - a*x1 + ax2rl * exp( r + a*ax2rl/b + b*x2 ) )
end

function tde_surflog(lx1, lx2, r, a, b)
    x1 = exp(lx1)
    x2 = exp(lx2)
    ax2rl = a*x2 - r + lx1 - lx2
    lx1 + r - a*x1 + ax2rl * exp( r + a*ax2rl/b + b*x2 )
end

n_star = 900
x_star = collect(range(0.001, length=n_star, stop=(13.0)))
y_star = collect(range(0.001, length=n_star, stop=(30.0)))
lx_star = collect(range(-2.0, length=n_star, stop=(3.25)))
ly_star = collect(range(-8.0, length=n_star, stop=(3.5)))

sys_mat = zeros(Float64, n_star, n_star)
logsys_mat = zeros(Float64, n_star, n_star)
tde_mat = zeros(Float64, n_star, n_star)
logtde_mat = zeros(Float64, n_star, n_star)
for i = 1:n_star
    for j = 1:n_star
        xx, tmp = sys_surf(x_star[i], y_star[j], r, a, b)
        global sys_mat[i,j] = isinf(xx) || abs(xx) > 15.0 ? NaN : xx

        xx, tmp = sys_surflog(lx_star[i], ly_star[j], r, a, b)
        global logsys_mat[i,j] = isinf(xx) || abs(xx) > 10.0 ? NaN : xx

        xx = tde_surf(x_star[i], x_star[j], r, a, b)
        global tde_mat[i,j] = isinf(xx) || abs(xx) > 15.0 ? NaN : xx

        xx = tde_surflog(lx_star[i], lx_star[j], r, a, b)
        global logtde_mat[i,j] = isinf(xx) || abs(xx) > 3.0 ? NaN : xx
    end
end

trace_sys = surface(Dict(
  :x => x_star,
  :y => y_star,
  :z => sys_mat,
  # :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))

trace_logsys = surface(Dict(
  :x => lx_star,
  :y => ly_star,
  :z => logsys_mat,
  # :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))

trace_tde = surface(Dict(
  :x => x_star,
  :y => x_star,
  :z => tde_mat,
  # :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))

trace_logtde = surface(Dict(
  :x => lx_star,
  :y => lx_star,
  :z => logtde_mat,
  # :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))

layout = Layout(Dict(
    :hovermode => "closest",  :margin => Dict( :r => 5, :t => 5, :b => 5, :l => 5),
    :scene => Dict(
        :xaxis => Dict(:title => "y[t-1]"),
        :yaxis => Dict(:title => "z[t-1]"),
        :zaxis => Dict(:title => "y[t]") )
    ))
data = [trace_sysdat, trace_sys]
plot(data, layout)

layout = Layout(Dict(
    :hovermode => "closest",
    :scene => Dict(
        :xaxis => Dict(:title => "log y[t-1]"),
        :yaxis => Dict(:title => "log z[t-1]"),
        :zaxis => Dict(:title => "log y[t]") )
    ))
data = [trace_logsysdat, trace_logsys]
plot(data, layout)

layout = Layout(Dict(
    :hovermode => "closest",
    :scene => Dict(
        :xaxis => Dict(:title => "y[t-1]"),
        :yaxis => Dict(:title => "y[t-2]"),
        :zaxis => Dict(:title => "y[t]") )
    ))
data = [trace_tdedat, trace_tde]
plot(data, layout)

layout = Layout(Dict(
    :hovermode => "closest",
    :scene => Dict(
        :xaxis => Dict(:title => "log y[t-1]"),
        :yaxis => Dict(:title => "log y[t-2]"),
        :zaxis => Dict(:title => "log y[t]") )
    ))
data = [trace_logtdedat, trace_logtde]
plot(data, layout)


text:
using Pkg
Pkg.activate("c:/git/JLBoost")
@time using DataFrames
@time using JDF
@time using JLBoost, LossFunctions
using JLBoost: LogitLogloss

###############################################################################
# fitting GiveMeSomeCredit
###############################################################################
if !isdir("c:/data/GiveMeSomeCredit/cs-training.jdf")
	using CSV
	a = CSV.read("c:/data/GiveMeSomeCredit/cs-training.csv", missingstrings=["NA"])
	rename!(a, Symbol("")=>:column1)
	create_missing!(a, :MonthlyIncome)
	create_missing!(a, :NumberOfDependents)
	type_compress!(a, compress_float=true)
	savejdf("c:/data/GiveMeSomeCredit/cs-training.jdf", a)
end

###############################################################################
# testing find_best_split
###############################################################################
@time a = loadjdf("c:/data/GiveMeSomeCredit/cs-training.jdf")

loss = LogitLogloss()

features = setdiff(names(a), [:SeriousDlqin2yrs])
target = :SeriousDlqin2yrs
warmstart = fill(0, length(a.age))
verbose = false
lambda = 0
gamma = 0

using ForwardDiff: derivative

f(warmstart, target) = begin
	p = 1/(1+exp(-(warmstart)))
	-(target*log(p) + (1-target)*log(1-p))
end

f(target) = warmstart -> f(warmstart, target)

x = derivative.(f.(a.SeriousDlqin2yrs), fill(0.0, nrow(a)))

x1 = deriv.(loss, a.SeriousDlqin2yrs, fill(0.0, nrow(a)))

x == x1


text:
# Julia wrapper for header: /usr/include/libavutil/timecode.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_timecode_adjust_ntsc_framenum2,
    av_timecode_get_smpte_from_framenum,
    av_timecode_make_string,
    av_timecode_make_smpte_tc_string,
    av_timecode_make_mpeg_tc_string,
    av_timecode_init,
    av_timecode_init_from_string,
    av_timecode_check_frame_rate

function av_timecode_adjust_ntsc_framenum2(framenum::Integer, fps::Integer)
    ccall((:av_timecode_adjust_ntsc_framenum2, libavutil), Cint, (Cint, Cint), framenum, fps)
end

function av_timecode_get_smpte_from_framenum(tc, framenum::Integer)
    ccall((:av_timecode_get_smpte_from_framenum, libavutil), UInt32, (Ptr{AVTimecode}, Cint), tc, framenum)
end

function av_timecode_make_string(tc, buf, framenum::Integer)
    ccall((:av_timecode_make_string, libavutil), Cstring, (Ptr{AVTimecode}, Cstring, Cint), tc, buf, framenum)
end

function av_timecode_make_smpte_tc_string(buf, tcsmpte::Integer, prevent_df::Integer)
    ccall((:av_timecode_make_smpte_tc_string, libavutil), Cstring, (Cstring, UInt32, Cint), buf, tcsmpte, prevent_df)
end

function av_timecode_make_mpeg_tc_string(buf, tc25bit::Integer)
    ccall((:av_timecode_make_mpeg_tc_string, libavutil), Cstring, (Cstring, UInt32), buf, tc25bit)
end

function av_timecode_init(tc, rate::AVRational, flags::Integer, frame_start::Integer, log_ctx)
    ccall((:av_timecode_init, libavutil), Cint, (Ptr{AVTimecode}, AVRational, Cint, Cint, Ptr{Cvoid}), tc, rate, flags, frame_start, log_ctx)
end

function av_timecode_init_from_string(tc, rate::AVRational, str, log_ctx)
    ccall((:av_timecode_init_from_string, libavutil), Cint, (Ptr{AVTimecode}, AVRational, Cstring, Ptr{Cvoid}), tc, rate, str, log_ctx)
end

function av_timecode_check_frame_rate(rate::AVRational)
    ccall((:av_timecode_check_frame_rate, libavutil), Cint, (AVRational,), rate)
end


text:
using Test

@testset "(32)-1->(22) fixed-direction signaling regression" begin

using BellScenario: Degeneracy, LocalPolytope, ConvexPolytope

# """
# maxwell_chitambar_inequalities()
#     Uses input/output relabelings to construct the complete set of bell inequalities
#     from a canonical subset of base inequalities.
#
# Outputs:
#     ineqs: Array, contains inequality vectors.
# """
function maxwell_chitambar_inequalities()

    joint_pos_ineq = [0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0]
    marginal_pos_ineq = [0 -1 0 0 1 0 0 0 0 0 0 0 0 0 0 0]
    I_ineq = [-1 -1 0 0 0 0 -1 1 1 0 -1 -1 0 0 1 0]
    II_ineq = [-1 0 0 0 -1 -1 -1 1 1 0 -1 -1 0 0 1 0]
    III_ineq = [-1 -1 -1 0 0 0 1 1 -1 1 -1 -1 0 0 1 -1]
    IV_ineq = [-1 0 0 0 -1 -1 -1 1 1 0 -1 -1 -1 1 1 0]
    V_ineq = [-1 -1 0 0 1 1 -1 -1 -1 1 0 0 -1 -1 1 -1]
    VI_ineq = [-1 -3 0 0 2 2 -2 1 0 -1 -1 -1 -1 0 1 -2]
    VII_ineq = [-1 -3 0 0 2 2 -2 1 0 -1 -1 -1 -2 1 1 -2]
    VIII_ineq = [-1 -3 0 0 2 2 -2 1 1 -2 -1 -1 -2 1 1 -2]

    ineqs = []

    # TODO: construct relabelngs in a way that is inclusive, but does not overshoot
    # by a factor of 4.5x duplicates.
    in_relabels = Degeneracy.bipartite_input_relabels((3,2),(2,2),"fixed-direction")
    (α_out_relabels, β_out_relabels) = Degeneracy.bipartite_output_relabels((3,2),(2,2),"fixed-direction")
    for in_relabel in in_relabels
        for α_out_relabel in α_out_relabels
            for β_out_relabel in β_out_relabels

                relabel = in_relabel*α_out_relabel*β_out_relabel

                push!(ineqs, joint_pos_ineq*relabel)
                push!(ineqs, marginal_pos_ineq*relabel)
                push!(ineqs, I_ineq*relabel)
                push!(ineqs, II_ineq*relabel)
                push!(ineqs, III_ineq*relabel)
                push!(ineqs, IV_ineq*relabel)
                push!(ineqs, V_ineq*relabel)
                push!(ineqs, VI_ineq*relabel)
                push!(ineqs, VII_ineq*relabel)
                push!(ineqs, VIII_ineq*relabel)
            end
        end
    end

    # TODO: fix so that we do not need to de-dupe
    unique(ineqs)
end

@testset "verifying polytope" begin

    dir = "./test/regression/files/"

    vertices = LocalPolytope.vertices((3,2),(2,2), bits=1)

    @test size(vertices) == (320,)
    @test size(vertices[1]) == (15,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)
    @test constraints[2] == []
    @test constraints[3] == []

    ineq_matches = maxwell_chitambar_inequalities()

    inequalities = constraints[1]
    @test issetequal(inequalities, ineq_matches)
end

end


text:
using DifferentialEquations, ReservoirComputing, Surrogates, ModelingToolkit, Plots, MLJLinearModels, SparseArrays
#=

	To generate a CTESN for a stiff ODE system:
		- specify system model ODE
		- specify a non-special initial condition
		- specify upper and lower bounds for all parameters in model
		- select set of parameter vectors sampled from parameter distribution, p[i] ~ U(p_lower[i], p_upper[i])
		- generate initial solutions using stiff ODE model for each p[i] => u(t,p[i])
		- generate ESN ODE
		- train ESN ODE on each stiff ODE solution, p[i] => u(t,p[i]) => WOut(p[i])
		- interpolate set of WOut(p[i]) using Surrogates
		- test CTESN at new q ~ in {p[i]}_i

=#

# define original stiff ODE system
function robertson!(du, u, p, t)
	du[1] = -p[1] * u[1] + p[3] * u[2] * u[3]
	du[2] =  p[1] * u[1] - p[3] * u[2] * u[3] - p[2] * u[2]^2
	du[3] =                                     p[2] * u[2]^2
	return nothing
end

# define original parameters and initial condition
p0 = (0.04, 3e7, 1e4)
u0 = [0.7; 0.2; 0.3]; u0 = u0./sum(u0) # to enforce sum(u) == 1.0
tspan = (0.0, 1e7)
tt = 10.0.^collect(range(-5.0, +5.0; length=33))
modelODESize = length(u0)

# define ODEProblem, optimize it, and solve for original ODE solution at (u0, p0)
prob = ODEProblem(robertson!, u0, tspan, p0)
#sol = solve(prob, Rosenbrock23(); abstol=1e-6, reltol=1e-6, saveat=tt)
sys = modelingtoolkitize(prob)
f   = eval(ModelingToolkit.generate_function(sys)[2])
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
prob = ODEProblem(ODEFunction(f, jac=jac), u0, tspan, p0)
sol = solve(prob, Rosenbrock23(); abstol=1e-6, reltol=1e-6, saveat=tt)

reservoirSize = 3000
iterative = reservoirSize > 300

W   = sparse(ReservoirComputing.init_reservoir_givendeg(reservoirSize, 0.8, 6))
Win = ReservoirComputing.init_dense_input_layer(reservoirSize, modelODESize, 0.1)

function CTESN!(du, u, p, t)
	du .= tanh.(p[1]*u .+ p[2]*p[3](t))
	return nothing
end
esnprob = ODEProblem(CTESN!, randn(Float64, reservoirSize), tspan, (W, Win, sol))
r = solve(esnprob; saveat = sol.t)

@assert isapprox(sol.t, r.t)

# compute first Wout
Wout = zeros(Float64, modelODESize, reservoirSize)
function fitData!(Wout, xData, rData; beta = 0.01, iterative = iterative)
	fitModel = RidgeRegression(beta)
	#=
	Note:
		rData = hcat(r(sol.t)...)
		xData = hcat(sol(sol.t)...)
	But we avoid re-interpolating them every time by requiring the r solution use `saveat=sol.t`
	=#
	for n=1:modelODESize
		tmp = transpose(fit(fitModel, transpose(rData[:,:]), xData[n,:]; solver=Analytical(iterative=iterative)))
		Wout[n,:] .= tmp[1:reservoirSize]
	end
	return nothing
end
fitData!(Wout, sol[:,:], r[:,:]; beta=0.01)

# define lower and upper bounds for each parameter
p_lower = [0.036, 2.7e7, 0.9e4]
p_upper = [0.044, 3.3e7, 1.1e4]

# select set of Sobol-sampled parameter vectors
p	= sample(100, p_lower, p_upper, SobolSample())
Wout 	= [zeros(Float64, modelODESize, reservoirSize) for _ in 1:length(p)]

# evaluate the ESN at each p, returning WOut(p)
function fitSamples!(Wout, p)

	function prob_func(prob, i, repeat)
		return remake(prob, p=p[i])
	end
	ensprob = EnsembleProblem(prob, prob_func = prob_func)
	sim = solve(ensprob, Rosenbrock23(), EnsembleThreads(); trajectories=length(p), abstol=1e-6, reltol=1e-6, saveat=sol.t)
	for n in 1:length(p)
		print("Fitting p[$(n)]...   \t")
		fitData!(Wout[n], sim[n], r; beta = 1e-6)
		print("Done.\n")
	end
	return nothing
end
fitSamples!(Wout, p)

# form interpolant over sampled parameter vectors
WoutInterpolant = RadialBasis(p, Wout, p_lower, p_upper)

# define unseen p̂ ensure it is not in p
p̂ = (0.0375,3.1e7, 1.01e4)
isinp = false;
for n in 1:length(p), m in 1:length(p̂)
	if isapprox(p̂[m], p[n][m])
		global isinp = true
		break
	end
end
@assert !isinp

# form Wout(p̂) and approximate x(t;p̂) from WoutInterpolant at unseen p̂
Ŵout = reshape(WoutInterpolant(p̂), size(Wout[1]))
x = Ŵout*r
#xInterpolant = RadialBasis(collect(transpose(sol.t)), x, [tspan[begin]], [tspan[end]])

sol_hat = solve(remake(prob, p=p̂), Rosenbrock23(); abstol=1e-6, reltol=1e-6, saveat=sol.t)
p1 = plot(sol)
plot!(p1, sol_hat, linestyle=:dash, linewidth=2)
plot!(p1, sol.t, transpose(x), linestyle=:dot, linewidth=3)
plot!(p1, xscale=:log)
p2 = plot(sol.t, transpose(sum(sol,dims=1)).-1.0)
plot!(p2, sol_hat.t, transpose(sum(sol_hat,dims=1)).-1.0, linestyle=:dash, linewidth=2)
plot!(p2, sol.t, transpose(sum(x,dims=1).-1.0), linestyle=:dot, linewidth=3)
plot!(p2, xscale=:log)
plot(p1,p2,layout = (2, 1))
savefig("./ODE_newODE_newCTESN.svg")

max_stars_repo_name: cmarcotte/CTESN

max_issues_repo_name: cmarcotte/CTESN

max_forks_repo_name: cmarcotte/CTESN


text:
using Makie, AbstractPlotting

function plot_solution(n::Int64,sol::Array{Float64,1},desiredstate::Function,alpha::Float64)
	# Simple scatter-plot of the nodes
	# Input: solution sol=(y,p,lambda)
	#        Desired state: y_d
	#	 L2-regularization alpha

	n_block::Int64 = Int64(length(sol)/3)
	n::Int64 = Int64( sqrt(n_block)-1 )

	dof_map = load("data"*string(n)*".jld2", "dof_map")

	xval::Array{Float64,1} = zeros(n_block)
	yval::Array{Float64,1} = zeros(n_block)

	stateval::Array{Float64,1} = zeros(n_block)
	controlval::Array{Float64,1} = zeros(n_block)
	adjointval::Array{Float64,1} = zeros(n_block)
	lagrangeval::Array{Float64,1} = zeros(n_block)
	y_d::Array{Float64,1} = zeros(n_block)

	# Separating all variables into vectors
	for i=1:n_block
		xval[i] = dof_map[i,:][1]
		yval[i] = dof_map[i,:][2]
		y_d[i] = desiredstate(xval[i],yval[i])
		stateval[i] = sol[i]
		adjointval[i] = sol[i+n_block]
		lagrangeval[i] = sol[i+2*n_block]
		controlval[i] = (1. /alpha)*(adjointval[i] - lagrangeval[i] )
	end

	# Plot state
	scene1 = Scene()
	scatter!(scene1 ,xval, yval, stateval, markersize=0.005, color = stateval ,colormap = :viridis)
	scale!(scene1, 1, 1, 1/(maximum(abs.(y_d))*2 ) )
	stateaxis  = scene1[Axis]
	stateaxis[:names, :axisnames] = ("x_1", "x_2", "y")

	# Plot Control
	scene2 = Scene()
	s = scatter!(scene2, xval, yval, controlval, markersize=0.005, color= controlval, colormap=:viridis )#color= :Black)
	scale!(scene2, 1, 1, 1/(maximum(abs.(controlval))*2 ) )
	controlaxis  = scene2[Axis]
	controlaxis[:names, :axisnames] = ("x_1", "x_2", "u")
	display(scene2)

	return scene1,scene2
end

max_stars_repo_name: dravinsi/demoPDEOPT

max_issues_repo_name: dravinsi/demoPDEOPT

max_forks_repo_name: dravinsi/demoPDEOPT