text:
################################
# Coded by Julien Pascal
# Last modification : 16/04/2016
# Replication of the the following paper:

####################################################################################
# "On the Dynamics of Unemployment and Wage Distributions"
# Paper from Jean Marc Robin (2011)
# Available here: https://sites.google.com/site/jmarcrobin/research/publications
####################################################################################

####################################################
# Model that performs the simulation of the economy
module simulation_economy

	using PyPlot, Distributions, DataFrames, StatsFuns, JLD, Copulas
	#to install the package "Copulas": Pkg.clone("https://github.com/floswald/Copulas.jl.git")
	#package from Florian Oswald

	export Simulate_Economy

	#################################
	# PATHS
	# change "path_main" if necessary
	path_main = "/home/julien/Final-Project-Julien-Pascal"

	path_table = string(path_main,"/tables/")  #the path tables
	path_surpluses = string(path_main,"/surpluses/")#path to precalculated objects, to speed up the process
	path_figures = string(path_main,"/figures/") #where to save the figures

	#Application from [0,1] intp [1,M]
	#To find the index corresponding to the percentile
	#
	#Use the fact that the Beta distribution is on [0,1]
	#and the function was discretized on an equidistant grid
	#
	# p = vector of pecentiles
	# M = integer, number of points
	function percentile_to_index(p::Float64,M::Int)

		index = 0 #initialization
		index = round(Int,(M-1)*p +1) #computation

		return index
	end

    ##################
    # Function that returns the inverse of lognormal distribution
    function logninv(p,mu,sigma)

        #output: scalar x
        logx0 = -sqrt(2)*erfcinv(2*p);
        x = exp(sigma*logx0 + mu);

        return x
    end
    #################

    ##########################
    #define match productivity:
    function yi_m(y_index,x_index, ygrid, xgrid)
        match = ygrid[y_index,1]*xgrid[x_index,1]
        return match
    end
    ##########################

    ##########################
    # define opportunity cost of employment:
    function zi_m(y_index,x_index, ygrid, xgrid, z0, alpha)

        op_cost = z0 + alpha*(yi_m(y_index,x_index, ygrid, xgrid) - z0) #call the function defined above

        return op_cost
    end

    #################################
    # Exit from unemployment rate: ft
    function ft_f(St_m,ut_m,ut,lm,lambda0,M)
    #St_m a vector of dimension M times 1

        ft = 0;
        for i=1:M
            #calculate the sum
            if St_m[1,i]> 0
                ft = ft + ut_m[1,i]*lm[i,1];
            end
        end
        #multiply by lamba0 and divide by the overall level of unemployment ut:
        ft = lambda0*ft/ut;
        return ft
    end

    ##################################################
    #function that calculates the job destruction rate
    function st_f(St_m,ut_m,ut,lm, delta, M)

        st = 0;
        for i=1:M
        #calculate the integral
            if St_m[1,i] <= 0
                st = st + (1 - ut_m[1,i])*lm[i,1];
            end
        end

        #finish the calculation:
        st = delta + (1 - delta)*st/(1 - ut);
        return st
    end

    #############################################################
    # function that calculate the measure of workers at hiw wages
    function gt_plus1_high_f(ut_m, gt_wi_low, gt_wi_high, i, lm, delta, lambda1, M, W_low, W_high, Si_m)

        gt_plus1 = zeros(1,M);

        for m = 1:M
            if Si_m[i,m]> 0
                #calculate the sum in the last paranthesis:
                sum = 0;
                if W_low[i,i,m] > Si_m[i,m]
                    sum = sum + gt_wi_low[1,m];
                end

                if W_high[i,i,m] > Si_m[i,m]
                    sum = sum + gt_wi_high[1,m];
                end

                gt_plus1[1,m] = (1 - delta)*(lambda1*(1 - ut_m[1, m])*lm[m,1] + (1 - lambda1)*((gt_wi_high[1,m] + sum)));
            end
        end
        return gt_plus1 # output = 1 times M vector
    end

    #############################################################
    #function that calculates the measure of workers at low wages
    function gt_plus1_low_f(ut_m, gt_wi_low, gt_wi_high, i, lm, lambda0, delta, lambda1, M, W_low, W_high, Si_m)
        #Inputs:
        # ut_m = 1 times M vector
        # gt_wi_low =  1 times M vector
        # gt_wi_high = 1 times M vector
        # i : integer, state of the economy

        gt_plus1 = zeros(1,M);

        for m =1:M
            #calculate the sum in the last paranthesis:
            if Si_m[i,m]> 0
                sum = 0;
                if W_low[i,i,m] < 0
                    sum = sum + gt_wi_low[1,m];
                end

                if W_high[i,i,m] < 0
                    sum = sum + gt_wi_high[1,m];
                end

                gt_plus1[1,m] = lambda0*ut_m[1,m]*lm[m,1] + (1 - delta)*(1 - lambda1)*(gt_wi_low[1,m] + sum);
            end
        end
        return gt_plus1 #Outuput is 1 times M vector
    end

    ##############################################
    #Calculates the next period's unemployment rate
    function ut_plus1_f(St_m, ut_m, delta, lambda0, M)

        #initialization:
        ut_plus1_m = zeros(1,M);

        for i=1:M
            #if the surplus is negative, unemployment is 1:
            if St_m[1,i] <= 0
                ut_plus1_m[1,i]= 1;
            else
                ut_plus1_m[1,i]= ut_m[1,i] + delta*(1 - ut_m[1,i]) - lambda0*ut_m[1,i];
            end
        end
        return ut_plus1_m
    end

    ########################################
    #function that calculates the high wages
    function wi_high_f(i, discount, Markov, lambda1, M, N, Si_m, W_high_star, zi_m_grid)

        #initialization
        wi_m_high = zeros(1,M);

        for m = 1:M
            #calculate the sum
            sum = 0;
            for j=1:N
                if Si_m[j,m]> 0
                sum = sum + Markov[i,j]*((lambda1*Si_m[j,m]) + (1 - lambda1)*W_high_star[j,i,m]);
                end
            end
            wi_m_high[1,m] = Si_m[i,m] + zi_m_grid[i,m]- discount*(sum);
        end
        return wi_m_high
    end

    ########################
    #Calculate the low wages:
    function wi_low_f(i, discount, Markov, lambda1, M, N, Si_m, W_low_star, zi_m_grid)

        #input: i = integer, state of the economy

        #initialization
        wi_m_low = zeros(1,M);
        for m = 1:M
            #calculate the sum
            sum = 0;
            for j=1:N
                if Si_m[j,m]> 0
                sum = sum + Markov[i,j]*((lambda1*Si_m[j,m]) + (1 - lambda1)*W_low_star[j,i,m]);
                end
            end
            wi_m_low[1,m] = zi_m_grid[i,m] - discount*(sum);
        end
        return wi_m_low
    end

    ############
    # Parameters
    ############

    #############################################
    # Values of the parameters:
    # A. Values from the 2011 paper
    # return a dictionary:
    function parameters_papers()
	    z0 = 0.77
	    sigma = 0.023
	    pho = 0.94
	    lambda0 = 0.99
	    k = 0.12
	    lambda1  = k*lambda0
	    s = 0.42
	    x_lower_bound = 0.73
	    eta = 2.00
	    mu = 5.56
	    alpha = 0.64
	    tau = 0.5
		delta = 0.042 # "4.2# exogenous layoff rate"
		r = 0.05/4 # quarterly interest rate
		discount = (1 - delta)/(1 + r)
		epsilon = 0.002
		N = 100 #number of states
		M = 500 #number of ability levels

	    return Dict("z0"=> z0, "sigma"=> sigma,"pho"=>pho,"lambda0"=>lambda0,"k"=>k, "lambda1"=>lambda1, "s"=>s, "x_lower_bound"=>x_lower_bound, "eta"=>eta, "mu"=>mu,"alpha"=>alpha, "tau" => tau, "delta" => delta, "r" =>r , "discount" =>discount, "epsilon" =>epsilon, "N"=>N, "M"=>M)
 	end

    #########################################
    # B. Estimated values for the parameters:
    # Read the coefficients estimated by the method of simulated moments:
    # Check the table exists
    function parameters_estimated()

	    # call it read estimated_parameters
	    if isfile(string(path_table, "estimated_parameters.csv"))
	    	println("loading estimated parameters")
			coeff = readtable(string(path_table, "estimated_parameters.csv"))
			"""
			z0 = coeff[:Value][1]
		    sigma = coeff[:Value][2]
		    pho = coeff[:Value][3]
		    lambda0 = coeff[:Value][4]
		    eta = coeff[:Value][5]
		    mu = coeff[:Value][6]
		    x_lower_bound = coeff[:Value][7]
		    delta = coeff[:Value][8] #exogenous layoff rate

		    k = 0.12
		    lambda1  = k*lambda0
		    #x_lower_bound = 0.73
		    alpha = 0.64
		    ##############################################""

		    tau = 0.5
		    #delta = 0.046
		    r = 0.05/4 # quarterly interest rate
		    discount = (1 - delta)/(1 + r)
		    epsilon = 0.002
		    """

		    """
		    z0 = coeff[:Value][1];
		    sigma = coeff[:Value][2];
		    lambda0 = coeff[:Value][3];
		    eta =  coeff[:Value][4];
		    mu =  coeff[:Value][5];
		    x_lower_bound = coeff[:Value][6] #lower bound of the ability level
		    alpha = coeff[:Value][7]
		    delta = coeff[:Value][8]

		    pho = 0.94
		    k = 0.12;
		    lambda1  = k*lambda0;
		    #x_lower_bound = 0.73
		    #alpha = 0.64
		    tau = 0.5
		    #delta = 0.043 # exogenous layoff rate

		    r = 0.05/4 # interest rate
		    discount = (1 - delta)/(1 + r)
		    epsilon = 0.002
		    """

			z0 = coeff[:Value][1];
		    sigma = coeff[:Value][2];
		    lambda0 = coeff[:Value][3];
		    eta =  coeff[:Value][4];
		    mu =  coeff[:Value][5];
		    x_lower_bound = coeff[:Value][6] #lower bound of the ability level

		    pho = 0.94
		    k = 0.12;
		    lambda1  = k*lambda0;
		    #x_lower_bound = 0.73
		    alpha = 0.64
		    tau = 0.5
		    delta = 0.042 # exogenous layoff rate

		    r = 0.05/4 # interest rate
		    discount = (1 - delta)/(1 + r)
		    epsilon = 0.002

		    N = 100 #number of states
			M = 500 #number of ability levels

			return Dict("z0"=> z0, "sigma"=> sigma,"pho"=>pho,"lambda0"=>lambda0,"k"=>k, "lambda1"=>lambda1, "x_lower_bound"=>x_lower_bound, "eta"=>eta, "mu"=>mu,"alpha"=>alpha, "tau" => tau, "delta" => delta, "r" =>r , "discount" =>discount, "epsilon" =>epsilon, "N"=>N, "M"=>M)

		else
			println(string(path_table, "estimated_parameters.csv"," not found"))
			println(string("parameters from the orginal paper were used instead"))
			output_dictionary = parameters_papers()
			return output_dictionary
		end

	end

	###########################################
	# function that generates the state space:
	# Return a dictionary with the different grids:
	# Elements: N, M, xgrid, agrid, ygrid, Markov, cop, b lm, zi_m_grid
	# Return the dictionary "parameters" as well
	function setup_grids()

		#Load parameters
		params = parameters_estimated()
		x_lower_bound = params["x_lower_bound"]
		epsilon = params["epsilon"]
		sigma = params["sigma"]
		pho= params["pho"]
		eta = params["eta"]
		mu = params["mu"]
		z0 = params["z0"]
		alpha = params["alpha"]
	    N = params["N"] #number of states
	    M = params["M"] #number of ability levels

	    # Define the grid for ability x
	    xgrid_space = linspace(x_lower_bound,x_lower_bound+1,M); # column vector with M rows

	    xgrid = Array{Float64}(M,1);
	    for i=1:M
	    	xgrid[i,1] = xgrid_space[i];
	    end

	    # Define the grid for match ability yi:
	    # 1. a:
	    agrid = linspace((0+epsilon),(1-epsilon),N);

	    # 2. yi:
	    #intialization:
	    ygrid = Array{Float64}(N,1);
	    for i=1:N
	    	ygrid[i] = logninv(agrid[i],0,sigma); #calculate the inverse of the lognormal
	    end

	    #Markov transition matrix:
	    Markov = zeros(N,N);

	    # number of dimension
	    # AR1 parameter
	    cop = Copulas.NormalCopula(2,pho)
	    #Gaussian copula pdf:
	    for j=1:N #move along the column:
	        for i=1:N #move along the lines:
	            Markov[i,j] = Copulas.dnormCopula([agrid[i] agrid[j]], cop)[];
	        end
	    end

	    #Normalize so that each row sum to 1:
	    for i = 1:N
	        Markov[i,:] = Markov[i,:]./sum(Markov[i,:]);
	    end

	    # distribution of workers on the grid x
	    lm = Array{Float64}(M,1);

	    # load the beta pdf
	    b = Beta(eta, mu)

	    for i=1:M
	        lm[i,1] = pdf(b, xgrid[i]-x_lower_bound); #evaluation of the beta pdf at a given point
	    end

	    #normalization so that mass of workers sum up to 1:
	    lm[:,1] = lm[:,1]./sum(lm[:,1]);

	    # Pre calculate the opportunity cost along the grid of ability:
	    # N times M matrix
	    # rows = aggregate staes
	    # columns = ability level
	    zi_m_grid = zeros(N,M);

	    for i=1:N
	        for m=1:M
	            zi_m_grid[i,m] = zi_m(i,m, ygrid, xgrid, z0, alpha);
	        end
	    end

		return Dict("xgrid"=>xgrid,"agrid"=>agrid,"ygrid"=>ygrid, "Markov"=>Markov, "cop"=>cop,
		"lm"=>lm, "b"=>b, "zi_m_grid"=>zi_m_grid, "params"=> params)
	end

    ###############################################
    # Load functions if already in the folder
    # Otherwise create them through value function iteration
    # Return a dictionary
    # Return the dictionaries "parameters" and "setup_grid" as well
    function create_surpluses_functions()

	    grids = setup_grids()
	    params = grids["params"]
	    N = params["N"] #number of states
	    M = params["M"] #number of ability levels
	    discount = params["discount"]

	    zi_m_grid = grids["zi_m_grid"]
	    xgrid = grids["xgrid"]
	    ygrid = grids["ygrid"]
	    Markov = grids["Markov"]

	    # Pre calculate values:
		G = zeros(N,M); #intialization
		p = zeros(N,M);
		for i = 1:N
			for m = 1:M
				G[i,m] = yi_m(i,m, ygrid, xgrid) - zi_m_grid[i,m];
				p[i,m] = yi_m(i,m, ygrid, xgrid);
			end
		end

    	# A. The match surplus
    	#=
	    if isfile(string(path_surpluses,"match_surplus.jld")) #load it if exists
	        println(string(path_surpluses,"match_surplus.jld"," found"))
	        println("loading it")
	        #Si_m = jldopen(string(path_surpluses,"match_surplus.jld"), "r") do file
	        #    read(file, "Si_m")
	        #end
	        file0 =  jldopen(string(path_surpluses,"match_surplus.jld"), "r")
		    Si_m = read(file0, "Si_m")
		    #close(file0)
	    else #otherwise create it by vale function iteration
	    =#

		    ##########################################################
		    # Calculate the match surplus by value function iteration:
		    # N times M matrices
		    # N = aggregate level
		    # M = ability level
		    println("Calculating the match surplus:")
	        Si_m = zeros(N,M);

	        #tol = 0.0001;
	        tol = 0.01;
	        maxits = 300;
	        dif = tol+tol;
	        its = 1;

	        #initialization:
	        up = ones(N,M);
	        up_plus1 = zeros(N,M);
	        compare = zeros(N,M);

		        while dif>tol
		            up_plus1 =  G + discount*Markov*max(up,compare)
		            dif = norm(up_plus1 - up)
		            up = up_plus1
		            its = its + 1
		            print(its)

		            if its > maxits
		                break
		            end
		        end

	        Si_m = up
	        #save the match surplus:
	       	#jldopen(string(path_surpluses,"/match_surplus.jld"), "w") do file
	        #	write(file, "Si_m", Si_m)
	        #end
	        #println(string(path_surpluses,"/match_surplus.jld")," created")
	    #end

	    # B. the wages:
		if isfile(string(path_surpluses,"W_low_star.jld"))
			println(string(path_surpluses,"W_low_star.jld", " found"))
			println("loading it")
		    #W_low_star = jldopen(string(path_surpluses,"W_low_star.jld"), "r") do file
		    #       read(file, "W_low_star")
		    #end

		    file1 =  jldopen(string(path_surpluses,"W_low_star.jld"), "r")
		    W_low_star = read(file1, "W_low_star")
		    #close(file1)

		else
			println(string(path_surpluses,"W_low_star.jld", " is missing"))
		    error("Run the function Wages() of the module calculate_wages")
		end

		if isfile(string(path_surpluses,"W_low.jld"))
			println(string(path_surpluses,"W_low.jld", " found"))
			println("loading it")
		    #W_low = jldopen(string(path_surpluses,"W_low.jld"), "r") do file
		    #       read(file, "W_low")
		    #end
		   	file2 =  jldopen(string(path_surpluses,"W_low.jld"), "r")
		    W_low = read(file2, "W_low")
		    #close(file2)

		else
			println(string(path_surpluses,"W_low.jld", " is missing"))
		    error("Run the function Wages() of the module calculate_wages")
		end

		if isfile(string(path_surpluses,"W_high_star.jld"))
			println("loading it")
			println(string(path_surpluses,"W_high_star.jld", " found"))
		    #W_high_star = jldopen(string(path_surpluses,"W_high_star.jld"), "r") do file
		    #        read(file, "W_high_star")
		    #end

		    file3 =  jldopen(string(path_surpluses,"W_high_star.jld"), "r")
		    W_high_star = read(file3, "W_high_star")
		    #close(file3)

		else
			println(string(path_surpluses,"W_high_star.jld", " is missing"))
		    error("Run the function Wages() of the module calculate_wages")
		end

		if isfile(string(path_surpluses,"W_high.jld"))
			println("loading it")
			println(string(path_surpluses,"W_high.jld", " found"))
		    #W_high = jldopen(string(path_surpluses,"W_high.jld"), "r") do file
		    #       read(file, "W_high")
		    #end

		    file4 =  jldopen(string(path_surpluses,"W_high.jld"), "r")
		    W_high = read(file4, "W_high")
		    #close(file4)

		else
			println(string(path_surpluses,"W_high.jld", " is missing"))
		    error("Run the function Wages() of the module calculate_wages")
		end

		return Dict("Si_m"=> Si_m ,"G"=>G,"p"=>p, "W_low_star"=>W_low_star, "W_low"=>W_low,
		"W_high_star"=>W_high_star, "W_high"=>W_high, "params" => params, "grids" => grids)
	end

	#############################################
	# Function that actually simulate the economy
	#############################################
	# returns a dictionary
	function execute_simulation(n_years = 2500)

		println("Simulate the economy")
		#Load the parameters, the grids and the value functions:
		Surpluses = create_surpluses_functions()

		params = Surpluses["params"]
		x_lower_bound = params["x_lower_bound"]
		epsilon = params["epsilon"]
		sigma = params["sigma"]
		pho= params["pho"]
		eta = params["eta"]
		mu = params["mu"]
		z0 = params["z0"]
		alpha = params["alpha"]
	    N = params["N"] #number of states
	    M = params["M"] #number of ability levels
	    lambda0 = params["lambda0"]
	    lambda1 = params["lambda1"]
	    delta = params["delta"]
	    discount = params["discount"]

	    grids = Surpluses["grids"]
	    zi_m_grid = grids["zi_m_grid"]
	    xgrid = grids["xgrid"]
	    ygrid = grids["ygrid"]
	    Markov = grids["Markov"]
	    lm = grids["lm"]
	    b = grids["b"] #beta distribution

	    Si_m = Surpluses["Si_m"]
	    G = Surpluses["G"]
	    p = Surpluses["p"]
	    W_low_star = Surpluses["W_low_star"]
	    W_low = Surpluses["W_low"]
	    W_high_star = Surpluses["W_high_star"]
	    W_high = Surpluses["W_high"]

	    # Load the randomly generated numbers:
		if isfile(string(path_surpluses,"random_numbers.jld")) #load it if exists
	        println(string(path_surpluses,"random_numbers.jld"," found"))
	        println("loading it")
	        file_random_numbers =  jldopen(string(path_surpluses,"random_numbers.jld"), "r")
		    random_numbers = read(file_random_numbers, "random_numbers")
	    else
	    	error("random numbers not generated")
	    end

		#n_years = 2500; # years to simulate
		number_periods = n_years*4; # one period = one quarter

		y_index_r = Array{Int16}(number_periods,1); #store the indexes of the shock
		y_r = Array{Float64}(number_periods,1) #store value of the shock

		#initial shock:
		y_index_r[1,1] = 50; #initial shock index
		y_r[1,1] = ygrid[y_index_r[1,1], 1]; #initial shock value

		#Initialization:
		ut_m_r = ones(number_periods,M); #row for the period, column for ability level
		St_m_r = Array{Float64}(number_periods,M); #row for the period, column for ability

		ut_r = Array{Float64}(number_periods,1);

		ft_r = Array{Float64}(number_periods,1);
		qt_r = Array{Float64}(number_periods,1);
		st_r = Array{Float64}(number_periods,1);

		wi_m_low_r = zeros(number_periods,M); #row for the period, column for ability;
		wi_m_high_r = zeros(number_periods,M); #row for the period, column for ability

		gt_wi_low_r = zeros(number_periods,M); #row for the period, column for ability;
		gt_wi_high_r = zeros(number_periods,M); #row for the period, column for ability;

		wi_m_r = zeros(number_periods,M);#row for the period, column for ability;
		wi_r = zeros(number_periods,1); #mean wage

		#measure of unemployed workers:
		m_unemployed_r = zeros(number_periods,M); #0 unemployed at t = 0

				#Indexes to find the percentiles in the distribution of workers:
		p_10 = percentile_to_index(quantile(b, 0.10), M) # index of the 10th percentile on the grid [1,M]
		p_20 = percentile_to_index(quantile(b, 0.20), M) # index of the 20th percentile on the grid [1,M]
		p_25 = percentile_to_index(quantile(b, 0.25), M) # etc.
		p_30 = percentile_to_index(quantile(b, 0.30), M)
		p_40 = percentile_to_index(quantile(b, 0.40 ),M)
		p_50 = percentile_to_index(quantile(b, 0.50), M)
		p_60 = percentile_to_index(quantile(b, 0.60), M)
		p_70 = percentile_to_index(quantile(b, 0.70), M)
		p_80 = percentile_to_index(quantile(b, 0.80), M)
		p_90 = percentile_to_index(quantile(b, 0.90), M)

		#unemployment by education:
		#25th percentile:
		u_25_p = zeros(number_periods,1);
		#50th percentile:
		u_50_p = zeros(number_periods,1);
		#90th percentile
		u_90_p = zeros(number_periods,1);

		#At t=1, half on the people at the starting salary and half of the people at the promotion salary
		gt_wi_low_r[1,:] = transpose(lm)/2
		gt_wi_high_r[1,:] = transpose(lm)/2

		########################
		# Loop over the economy:
		for t=1:(number_periods-1)

		    #println(string("quarter #:", t))
		    # Calculate the aggregate unemployment
		    ut_r[t,1] = dot(ut_m_r[t,:],lm[:,1]); #dot product.

		    # Calculate the surplus given the actual value of the shock
		    St_m_r[t,:] = Si_m[y_index_r[t,1],:];

		    #Remark: everytime I select a row, Julia treats it as a column vector
		    #That is why I have to transpose my inputs in the next function
		    # Exit rate from unemployment:
		    ft_r[t,1] = ft_f(St_m_r[t,:]', ut_m_r[t,:]', ut_r[t,1], lm, lambda0, M);

		    #Job destruction rate:
		    st_r[t,1] = st_f(St_m_r[t,:]', ut_m_r[t,:]', ut_r[t,1], lm, delta, M);

		    # Law of motion of unemployment:
		    ut_m_r[t+1,:] = ut_plus1_f(St_m_r[t,:]', ut_m_r[t,:]', delta, lambda0, M);

		    #Calculate group unemployment rate:
		    #25th percentile:
		    u_25_p[t+1,1] = dot(lm[1:p_25,1], ut_m_r[t+1,1:p_25])/sum(lm[1:p_25,1]);

		    #50th percentile:
		    u_50_p[t+1,1] = dot(lm[1:p_50,1], ut_m_r[t+1,1:p_50])/sum(lm[1:p_50,1]);

		    #90th percentile
		    u_90_p[t+1,1]= dot(lm[1:p_90,1], ut_m_r[t+1,1:p_90])/sum(lm[1:p_90,1]);

		    #Calculate the wages:
		    # 1. starting wages
		    wi_m_low_r[t,:] = wi_low_f(y_index_r[t,1], discount, Markov, lambda1, M, N, Si_m, W_low_star, zi_m_grid); # input: the index state of the economy

		    #2. promotion wages
		    wi_m_high_r[t,:] = wi_high_f(y_index_r[t,1], discount, Markov, lambda1, M, N, Si_m, W_high_star, zi_m_grid); #input: the index state of the economy

		    #measure of workers of ability m employed at low wage at the end of
		    #period t
		    gt_wi_low_r[t+1,:] = gt_plus1_low_f(ut_m_r[t,:]',gt_wi_low_r[t,:]',gt_wi_high_r[t,:]',y_index_r[t,1], lm, lambda0, delta, lambda1, M, W_low, W_high, Si_m);

		    #measure of workers of ability m employed at high wage at the end of
		    #period t
		    gt_wi_high_r[t+1,:] = gt_plus1_high_f(ut_m_r[t,:]', gt_wi_low_r[t,:]', gt_wi_high_r[t,:]', y_index_r[t,1], lm, delta, lambda1, M, W_low, W_high, Si_m);

		    #measure of unemployed workers at the end of period t:
		    # = mesure of workers with ability m minus employed people
		    m_unemployed_r[t+1,:] = (transpose(lm[:,1]) - (gt_wi_low_r[t+1,:]' + gt_wi_high_r[t+1,:]'));

		    #calculate the average wage by worker type:
		    #weight by the measure of workers with low/high wages:
		    for i=1:M #loop over ability level
		    	mass = gt_wi_low_r[t+1,i] + gt_wi_high_r[t+1,i]
		    	if (mass!= 0) & (isnan(mass)==false) #cannot divide by zero and don't want to deal with NaN
					wi_m_r[t,i] = (wi_m_low_r[t,i]*gt_wi_low_r[t+1,i] + wi_m_high_r[t,i]*gt_wi_high_r[t+1,i])/(gt_wi_low_r[t+1,i] + gt_wi_high_r[t+1,i])
				else
			   		wi_m_r[t,i] = NaN
			   	end
			end
		    ############################
		    # mean wage denoted by wi_r
		    #
		    # calculate recursively the mean wage:
			wi_r[t,1]= 0; #initialization
			csum = 0; #initialisation

			for i=1:M #loop over ability
				mass = gt_wi_low_r[t+1,i] + gt_wi_high_r[t+1,i]
				if (mass!= 0) & (isnan(mass)==false) & (isnan(wi_m_r[t,i])==false) #cannot divide by zero and don't want to deal with NaN
					wi_r[t,1] = wi_r[t,1] + wi_m_r[t,i]*mass; #weight by the measure of workers at the given wage
					csum = csum + mass;
				end
			end

			if csum != 0
				wi_r[t,1] = wi_r[t,1]/csum; #divide by the mass of workers who have a wage
			else
				wi_r[t,1] = NaN #no one has a wage.
			end
		    ######################

		    # New shock from the markov transition matrix:
		    #r = rand(); #draw a number from the uniform distribution on 0, 1
		    r = random_numbers[t]

		    # I use the Markov transition matrix previously defined:
		    prob = Markov[y_index_r[t,1],:];

		    #I stock the index in the Markov matrix, as well as the value of the
		    #shock y
		    y_index_r[t+1,1] = sum(r .>= cumsum([0; prob]));

		    y_r[t+1,1] = ygrid[y_index_r[t+1,1], 1];
		end

		#Surpluses: contains the surplus value functions, the parameters, the parameter grids
		#n_years: number years the simulation runs
		#number_periods: number of periods the simulation runs (one period = one quarter)
		#y_r: deviation of productivity from its long time tredn

		return Dict("Surpluses"=> Surpluses ,"n_years"=>n_years, "number_periods"=>number_periods, "y_r"=>y_r , "ut_m_r"=>ut_m_r,
		"ut_r"=>ut_r, "ft_r"=>ft_r, "st_r" => st_r, "wi_r" => wi_r, "wi_m_r" =>wi_m_r, "u_25_p"=>u_25_p, "u_50_p"=>u_50_p, "u_90_p"=>u_90_p,
		"gt_wi_low_r"=>gt_wi_low_r, "gt_wi_high_r"=>gt_wi_high_r, "lm"=>lm)
	end

	############################
	# Analysis of the simulation
	############################
	# model: has to be the output of the function execute_simulation()
	# get_rid_of: to discard the impact of initial condition, get rid of a given percentage of first observations
	function analyse_economy(model, get_rid_of = 10)

		println("Analyse the economy")
		Surpluses = model["Surpluses"]

		params = Surpluses["params"]
		x_lower_bound = params["x_lower_bound"]
		epsilon = params["epsilon"]
		sigma = params["sigma"]
		pho= params["pho"]
		eta = params["eta"]
		mu = params["mu"]
		z0 = params["z0"]
		alpha = params["alpha"]
	    N = params["N"] #number of states
	    M = params["M"] #number of ability levels
	    lambda0 = params["lambda0"]
	    lambda1 = params["lambda1"]
	    delta = params["delta"]
	    discount = params["discount"]

	    grids = Surpluses["grids"]
	    zi_m_grid = grids["zi_m_grid"]
	    xgrid = grids["xgrid"]
	    ygrid = grids["ygrid"]
	    Markov = grids["Markov"]
	    lm = grids["lm"]

	    Si_m = Surpluses["Si_m"]
	    G = Surpluses["G"]
	    p = Surpluses["p"]
	    W_low_star = Surpluses["W_low_star"]
	    W_low = Surpluses["W_low"]
	    W_high_star = Surpluses["W_high_star"]
	    W_high = Surpluses["W_high"]

	    n_years = model["n_years"]
	   	number_periods = model["number_periods"]
	   	y_r =  model["y_r"]
	   	ut_m_r = model["ut_m_r"]
	   	ut_r = model["ut_r"]
	   	ft_r = model["ft_r"]
	   	st_r = model["st_r"]
	   	wi_r = model["wi_r"]
	   	wi_m_r =  model["wi_m_r"]

	   	u_25_p = model["u_25_p"]
	   	u_50_p = model["u_50_p"]
	   	u_90_p = model["u_90_p"]

	   	gt_wi_low_r = model["gt_wi_low_r"]
	   	gt_wi_high_r =model["gt_wi_high_r"]

		discard = round(Int,floor(number_periods/get_rid_of)); #get rid of the first "get_rid_of"th observations

		#disrtribution of workers
		figure("Distribution_workers 1",figsize=(10,10))
		plot(xgrid[:],lm[:])
		xlabel("ability level")
		ylabel("mass of worker at each ability level")
		savefig("figures/Distribution_workers_by_ability.png")

		figure("Distribution_workers 2",figsize=(10,10))
		plot(lm[:])
		xlabel("Index ability level on [1,M]")
		ylabel("mass of worker at each index")

		# wage per decile:
		    # rows = time period
		    # columns = percentile. ex: first row = 10th decile, 2nd row= 20th decile
		w_p = zeros(number_periods,9); #initialization

		#Plot the surplus function:
        xgrid_plot = repmat(xgrid',N,1)
        ygrid_plot = repmat(ygrid,1,M)

        """
        fig = figure("pyplot_surfaceplot",figsize=(10,10))
        ax = fig[:add_subplot](2,1,1, projection = "3d")
        ax[:plot_surface](xgrid_plot, ygrid_plot, Si_m, rstride=2, edgecolors="k", cstride=2, cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
        xlabel("X")
        ylabel("Y")
        title("Surface Plot")

        subplot(212)
        ax = fig[:add_subplot](2,1,2)
        cp = ax[:contour](xgrid_plot, ygrid_plot, Si_m, colors="black", linewidth=2.0)
        ax[:clabel](cp, inline=1, fontsize=10)
        xlabel("X")
        ylabel("Y")
        title("Contour Plot")
        tight_layout()
        savefig("figures/Surplus_function.png")
        """

        fig = figure("Surplus1",figsize=(10,10))
        ax = fig[:add_subplot](1,1,1)
        cp = ax[:contour](xgrid_plot, ygrid_plot, Si_m, linewidth=2.0)
        ax[:clabel](cp, inline=1, fontsize=10)
        xlabel("ability xm")
        ylabel("productivity y")
        title("Contour Plot")
        tight_layout()
        savefig("figures/Surplus_function_1.png")

        fig = figure("Surplus2",figsize=(10,10))
        ax = fig[:add_subplot](1,1,1, projection = "3d")
        ax[:plot_surface](xgrid_plot, ygrid_plot, Si_m, rstride=2, edgecolors="k", cstride=2, cmap=ColorMap("gray"), alpha=0.8, linewidth=0.25)
        xlabel("X")
        ylabel("Y")
        title("Surface Plot")
         savefig("figures/Surplus_function_2.png")

		# Plot productivity match:
		figure("1",figsize=(10,10))
		for m=1:10:M
		    plot(ygrid[:],p[:,m])
		end
		xlabel("Aggregate State")
		ylabel("Productivity")
		#title("Match Productivity by Ability")
		savefig("figures/Simulation_Match_Productivity_by_Ability.png")

		figure("2",figsize=(10,10))
		time_plot = discard:(number_periods-1);
		# Plot the shock:
		plot(time_plot[:] ,y_r[discard:(number_periods-1),1])
		xlabel("Periods")
		ylabel("Aggregate shock")
		title("Aggregate shock over time")
		savefig("figures/Simulation_Aggregate_shock_over_time.png")

		figure("3",figsize=(10,10))
		#Plot aggregate unemployment:
		plot(time_plot[:],ut_r[discard:(number_periods-1),1])
		xlabel("Periods")
		ylabel("Unemployment rate")
		title("Aggregate Unemployment rate over time")
		savefig("figures/Simulation_Aggregate_Unemployment_rate_over_time.png")

		figure("4",figsize=(10,10))
		#Plot exit rate from unemployment:
		plot(time_plot[:],ft_r[discard:(number_periods-1),1])
		xlabel("Periods")
		ylabel("ft")
		title("Exit rate from unemployment")
		savefig("figures/Simulation_Exit_rate_from_unemployment.png")

		figure("5",figsize=(10,10))
		#Plot job destruction rate:
		plot(time_plot[:], st_r[discard:(number_periods-1),1])
		xlabel("Periods")
		ylabel("st")
		title("Job destruction rate")
		savefig("figures/Simulation_Job_destruction_rate.png")

		#plot the average wage
		figure("6",figsize=(10,10))
		plot(time_plot[:], wi_r[discard:(number_periods-1),1]);
		title("Mean Wage")
		xlabel("Periods")
		ylabel("wage")
		savefig("figures/Simulation_Mean_wage.png")

		#Unemployment rate by education level
		figure("7",figsize=(10,10))
		# 25% low skills
		scatter(ut_r[discard:(number_periods-1),1], u_25_p[discard:(number_periods-1),1], alpha=0.4, color="green", label ="25th percentile ability")
		# 50% low skills
		scatter(ut_r[discard:(number_periods-1),1], u_50_p[discard:(number_periods-1),1], alpha=0.7, color="red", label ="50th percentile ability")
		# 90% low skills
		scatter(ut_r[discard:(number_periods-1),1], u_90_p[discard:(number_periods-1),1], alpha=1, label ="90th percentile ability")
		title("Unemployment Rate Among Various Skill Groups")
		xlabel("overall unemployment rate")
		ylabel("unemployment rate by ability")
		legend(loc="best",fancybox="true")
		savefig("figures/Simulation_unemployment_rate_by_skill.png")

		#Scatterplot aggregate unemployment rate vs productivity
		figure("7bis",figsize=(10,10))
		scatter(y_r[discard:(number_periods-1),1], ut_r[discard:(number_periods-1),1], alpha=0.4)
		title("Productivity and Unemployment Rate")
		xlabel("Productivity")
		ylabel("Aggregate unemployment rate")
		legend(loc="best",fancybox="true")
		savefig("figures/Simulation_Productivity_and_Unemployment_Rate.png")

		"""
		# This part is quite slow (about 10 minutes)
		# Uncomment to get the dynamics of wage deciles
		#########################
		#Dynamics of wage decile:
		#########################

		###############################################
		# Calculate wage deciles and interdecile ratios

		trim = floor(Int16,1*discard); #to gain time, calculate for fewer periods

		D5_D1 = zeros(number_periods,1);
		D9_D1 = zeros(number_periods,1);
		D9_D5 = zeros(number_periods,1);

		for t = (number_periods-trim):(number_periods-1)
		println(t)

		    per = 0; #intialization

		    for i=1:9

		        #choose the decile of interest
		        per = per + 0.1; #10th, 20th, etc...
		        a = 0; #intialization
		        sump = 0;

		        while ((sump<per) & (a<500)) #find the ability index that corresponds to the decile "per"

					#initialize values:
					a = a+1;
					sump=0;
					csum=0;

		            #calculate a sum
		            for m = 1:a #loop over ability levels
						if ((wi_m_r[t,m] < wi_m_r[t,a]) & (wi_m_r[t,m]!=0) & (isnan(wi_m_r[t,m]) == false)) #take into account only employed workers
						sump = sump + (gt_wi_low_r[t+1, m] + gt_wi_high_r[t+1, m]); #weight by the distribution of workers
						end
		             end
					csum = sum(gt_wi_low_r[t+1,:] + gt_wi_high_r[t+1,:]);

					if csum !=0
						sump = sump/csum;
					else
						sump = 0;
					end
		        end

		        #store the value of the wage decile:
		        w_p[t,i] = wi_m_r[t,a];
		    end

		    #calculate interdecile ratio
		    D5_D1[t,1] = w_p[t,5]/w_p[t,1];
		    D9_D1[t,1] = w_p[t,9]/w_p[t,1];
		    D9_D5[t,1] = w_p[t,9]/w_p[t,5];
		end

		time_plot2 = (number_periods-trim):(number_periods-1);
		#Plot wage deciles:
		figure("8",figsize=(10,10))
		for m = 1:9 #loop over the deciles
		    plot(time_plot2[:], w_p[(number_periods-trim):(number_periods-1),m], label = string(m,"0th percentile"))
		end
		xlabel("Periods")
		ylabel("wage")
		title("Dynamics of Wage Deciles")
		legend(loc="best",fancybox="true")
		savefig("figures/Simulation_Dynamics_Wage_Deciles.png")

		#Plot Interdecile ratios:
		figure("9",figsize=(10,10))
		plot(time_plot2[:], D5_D1[(number_periods-trim):(number_periods-1),1], color="green", label = "D5_D1")
		plot(time_plot2[:], D9_D1[(number_periods-trim):(number_periods-1),1], color="red", label = "D9_D1")
		plot(time_plot2[:], D9_D5[(number_periods-trim):(number_periods-1),1], label = "D9_D5")
		xlabel("Periods")
		ylabel("Inter-decile ratios")
		title("Wage Inequalities")
		legend(loc="best",fancybox="true")
		savefig("figures/Simulation_Wage_Inequalities.png")

		#Apply the same transformations as for the empirical data:
		#Deviation from the mean (there is not trend in this setting, no need to claculate OLS)
		mean_wage = zeros(1,9); #column for percentile: first column = 10th percentile, etc... , 9th column = 90th percentile
		length_deviation = (number_periods-1) - (number_periods-trim) +1
		deviation_mean = zeros(length_deviation, 9) #row for time, column for ability;
		for m = 1:9
			#calculation:
			mean_wage[1, m] = mean(w_p[(number_periods-trim):(number_periods-1), m])
			deviation_mean[1:length_deviation, m] =  w_p[(number_periods-trim):(number_periods-1), m] - ones(length_deviation, 1)*mean_wage[1, m]
		end

		#Plot deviation from mean:
		figure("10",figsize=(10,10))
		for m = 1:9 #loop over the deciles
		    plot(time_plot2[:], deviation_mean[1:length_deviation, m],label = string(m,"0th percentile"))
		end

		xlabel("Periods")
		ylabel("wage")
		title("Deviations of Wages Deciles from Mean")
		legend(loc="best",fancybox="true")
		savefig("figures/Simulation_Deviation_Wage_Deciles_From_Mean.png")

		#Calculate the standard deviation of wage deciles:
		std_wage_deciles = zeros(1, 9) #column for percentile: first column = 10th percentile, etc... , 9th column = 90th percentile

		for m = 1:9 #loop over the percentiles
		    std_wage_deciles[1, m] = std(deviation_mean[1:length_deviation, m])
		end

		params_wages = ["10th", "20th", "30th", "40th", "50th", "60th", "70th", "80th", "90th"]
		wages_df = DataFrame()
		wages_df[:Parameters] = params_wages

		println(transpose(round(std_wage_deciles,5)))

		wages_df[:Stdev_Wage_Decile] = 0.0 #initializaion

		for m =1:9
			wages_df[m,:Stdev_Wage_Decile] = sum(round(std_wage_deciles[1, m],5))
		end

		println("Simulation results, Wages:")
		println(wages_df)
		writetable("tables/results_simulation_wages.csv", wages_df)

		#Summary statistics:
		mean_simulation = zeros(5)
		mean_simulation[1] = sum(mean(y_r[discard:(number_periods-1),1]))
		mean_simulation[2] = sum(mean(ut_r[discard:(number_periods-1),1]))
		mean_simulation[3] = sum(mean(ft_r[discard:(number_periods-1),1]))
		mean_simulation[4] = sum(mean(st_r[discard:(number_periods-1),1]))
		mean_simulation[5] = sum(mean(wi_r[discard:(number_periods-1),1]))

		std_simulation = zeros(5)
		std_simulation[1] = sum(std(log(y_r[discard:(number_periods-1),1])))
		std_simulation[2] = sum(std(log(ut_r[discard:(number_periods-1),1])))
		std_simulation[3] = sum(std(log(ft_r[discard:(number_periods-1),1])))
		std_simulation[4] = sum(std(log(st_r[discard:(number_periods-1),1])))
		std_simulation[5] = sum(std(log(wi_r[discard:(number_periods-1),1])))

		skewness_simulation = zeros(5)
		skewness_simulation[1] = sum(skewness(log(y_r[discard:(number_periods-1),1])))
		skewness_simulation[2] = sum(skewness(log(ut_r[discard:(number_periods-1),1])))
		skewness_simulation[3] = sum(skewness(log(ft_r[discard:(number_periods-1),1])))
		skewness_simulation[4] = sum(skewness(log(st_r[discard:(number_periods-1),1])))
		skewness_simulation[5] = sum(skewness(log(wi_r[discard:(number_periods-1),1])))

		kurtosis_simulation = zeros(5)
		kurtosis_simulation[1] = sum(kurtosis(log(y_r[discard:(number_periods-1),1])))
		kurtosis_simulation[2] = sum(kurtosis(log(ut_r[discard:(number_periods-1),1])))
		kurtosis_simulation[3] = sum(kurtosis(log(ft_r[discard:(number_periods-1),1])))
		kurtosis_simulation[4] = sum(kurtosis(log(st_r[discard:(number_periods-1),1])))
		kurtosis_simulation[5] = sum(kurtosis(log(wi_r[discard:(number_periods-1),1])))

		params_name = ["y", "u", "ft", "st", "w"]
		df = DataFrame()
		df[:Parameters] = params_name
		df[:Mean] = round(mean_simulation, 4)
		df[:Std] = round(std_simulation, 4)
		df[:Skewness] = round(skewness_simulation, 4)
		df[:Kurtosis] = round(kurtosis_simulation, 4)

		println("Simulation results, Moments:")
		println(df)
		writetable("tables/results_simulation_moments.csv", df)

		return Dict("Moments" => df, "Wages_deciles"=> wages_df)
		"""

	end

	#Function that simulates the economy based on
	#the estimated parameters:
	function Simulate_and_Anlyse_Economy()

		#exectue the simulation:
		model_ouput = execute_simulation()
		analyse_economy(model_ouput, 10) #get rid of the first 10th observations

		println("enter y to close this session.")
		ok = readline(STDIN) # Equivalent of the function "input" on the version of Julia I am running
		if ok == "y"
			quit()
		end
	end
end


text:
@testset "geometric sampler" begin

	cc = dummyChoiceContext()
	cc.derivationstate.currentrecursiondepth = 3

	@testset "getparams, numparams, paramranges, sample" begin

		@testset "using default construction" begin

		    s = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3)

		    @testset "numparams and paramranges" begin
		        @test DataGenerators.numparams(s) == 3
		        prs = DataGenerators.paramranges(s)
		        @test typeof(prs) <: Vector{Tuple{Float64,Float64}}
		        @test prs == [(0.0,1.0),(0.0,1.0),(0.0,1.0)]
		    end

		    @testset "default params" begin
		        @test DataGenerators.getparams(s) == [0.5, 0.5, 0.5]
		    end

			@mtestset "default sampling consistent with base sampler depth 3" reps=Main.REPS alpha=Main.ALPHA begin
				x, trace = DataGenerators.sample(s, (0,1), cc)
				@test typeof(x) <: Int
				@mtest_values_include [0,1,2] x
				@mtest_distributed_as Geometric(0.5) x
			end

		end

		@testset "using non-default construction" begin

		    s = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3, [0.5, 0.1, 0.7])
		    @test DataGenerators.getparams(s) == [0.5, 0.1, 0.7]

		    @mtestset "default sampling consistent with base sampler depth 3" reps=Main.REPS alpha=Main.ALPHA begin
		    	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Geometric(0.7) x
		    end

		end

	end

	@testset "getparams, numparams, paramranges, sample using multi-parameter base sampler" begin

		@testset "using default construction" begin

		    s = DataGenerators.RecursionDepthSampler(DataGenerators.CategoricalSampler(4), 2)

		    @testset "numparams and paramranges" begin
		        @test DataGenerators.numparams(s) == 8
		        prs = DataGenerators.paramranges(s)
		        @test typeof(prs) <: Vector{Tuple{Float64,Float64}}
		        @test prs == [(0.0,1.0),(0.0,1.0),(0.0,1.0),(0.0,1.0),(0.0,1.0),(0.0,1.0),(0.0,1.0),(0.0,1.0),]
		    end

		    @testset "default params" begin
		        @test DataGenerators.getparams(s) == [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25,]
		    end

			@mtestset "consistent with base sampler max depth" reps=Main.REPS alpha=Main.ALPHA begin
				x, trace = DataGenerators.sample(s, (0,1), cc)
				@test typeof(x) <: Int
				@mtest_values_are [1,2,3,4] x
				@mtest_distributed_as Categorical(4) x
			end

		end

		@testset "using non-default construction" begin

		    s = DataGenerators.RecursionDepthSampler(DataGenerators.CategoricalSampler(4), 2, [0.25, 0.25, 0.25, 0.25, 0.1, 0.2, 0.3, 0.4,])
		    @test DataGenerators.getparams(s) == [0.25, 0.25, 0.25, 0.25, 0.1, 0.2, 0.3, 0.4,]

		    @mtestset "consistent with base sampler max depth" reps=Main.REPS alpha=Main.ALPHA begin
		    	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Categorical([0.1, 0.2, 0.3, 0.4,]) x
		    end

		end

	end

	@testset "setparams" begin

		s = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3)

	    @testset "valid parameters $params" for params in [[0.8, 0.4, 0.8,], [0.5, 0.1, 0.5], [0.7, 0.4, 0.6],]

	        DataGenerators.setparams!(s, params)

			@mtestset "consistent with base sampler depth 3" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Geometric(params[3]) x
			end

	    end

	    @testset "wrong number of parameters" begin

	        @test_throws ErrorException DataGenerators.setparams!(s, [0.8, 0.4,])
	        @test_throws ErrorException DataGenerators.setparams!(s, [0.8, 0.4, 0.1, 0.3])

	    end

	end

	@testset "setparams using multi-parameter base sampler" begin

		s = DataGenerators.RecursionDepthSampler(DataGenerators.CategoricalSampler(2), 3)

	    @testset "valid parameters $params" for params in [[0.5, 0.5, 0.5, 0.5, 0.5, 0.5,], [0.5, 0.5, 0.2, 0.8, 0.4, 0.6,], [0.3, 0.7, 0.1, 0.9, 0.7, 0.3,],]

	        DataGenerators.setparams!(s, params)

			@mtestset "consistent with base sampler depth 3" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Categorical(params[5:6]) x
			end

	    end

	    @testset "wrong number of parameters" begin

	        @test_throws ErrorException DataGenerators.setparams!(s, [0.5, 0.5, 0.5, 0.5, 0.5,])
	        @test_throws ErrorException DataGenerators.setparams!(s, [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,])

	    end

	end

	@testset "estimateparams" begin

		params = [0.7, 0.4, 0.7,]
	    otherparams = [0.5, 0.7, 0.5,]

	    @testset "from sufficient traces" begin

			s1 = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3, params)
			cc1 = dummyChoiceContext()

			alltraces = Any[]
			for depth in 1:3
				cc1.derivationstate.currentrecursiondepth = depth
		        traces = map(1:100) do i
		            x, trace = DataGenerators.sample(s1, (0,1), cc1)
		            trace
		        end
				alltraces = [alltraces; traces]
			end

			s2 = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3, otherparams)
		    estimateparams!(s2, alltraces)

			cc2 = dummyChoiceContext()

			@mtestset "consistent with otherparams for depth $depth" reps=Main.REPS alpha=Main.ALPHA for depth in 1:3
				cc2.derivationstate.currentrecursiondepth = depth
	        	x, trace = DataGenerators.sample(s2, (0,1), cc2)
		        @mtest_distributed_as Geometric(params[depth]) x
			end

	    end

	    @testset "too few traces" begin

			s1 = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3, params)
			cc1 = dummyChoiceContext()

			alltraces = Any[]
			for depth in [1,3]
				cc1.derivationstate.currentrecursiondepth = depth
		        traces = map(1:100) do i
		            x, trace = DataGenerators.sample(s1, (0,1), cc1)
		            trace
		        end
				alltraces = [alltraces; traces]
			end

			s2 = DataGenerators.RecursionDepthSampler(DataGenerators.GeometricSampler(), 3, otherparams)
		    estimateparams!(s2, alltraces)

			cc2 = dummyChoiceContext()

			@mtestset "consistent with otherparams for depth $depth" reps=Main.REPS alpha=Main.ALPHA for depth in 2
				cc2.derivationstate.currentrecursiondepth = depth
	        	x, trace = DataGenerators.sample(s2, (0,1), cc2)
		        @mtest_distributed_as Geometric(otherparams[depth]) x
			end

	    end

	end

end

ext: jl


text:
# See https://en.wikipedia.org/wiki/Padovan_sequence
# Inspired by https://math.stackexchange.com/questions/3798107/is-there-a-better-way-to-calculate-padovan-sequence-than-on

padovan(n) = ([0 1 0; 0 0 1; 1 1 0]^n)[3, 2]

expressions = [
    :(padovan(1))
    :(padovan(2))
    :(padovan(3))
    :(padovan(4))
    :(padovan(5))
    :(padovan(6))
    :(padovan(7))
    :(padovan(8))
    :(padovan(9))
    :(padovan(10))
    :(padovan(91))
    :(padovan(big(283)))
    :(padovan(big(3674)))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
####
# UM Oracle tests
####

facts("suppFcnTest UM") do
	lquants = [-.05, -.25]
	uquants = [.1, .2]

	w = UMOracle(lquants, uquants, 1e-6, false)
	zstar, ustar = suppFcn([1, 1], w, :Min)
	@fact zstar -->  roughly(-0.3)
	@fact ustar[1] --> roughly(-0.05)
	@fact ustar[2] --> roughly(-0.25)

	zstar, ustar = suppFcn([1, 1], w, :Max)
	@fact zstar --> roughly(.3)
	@fact ustar[1] --> roughly(.1)
	@fact ustar[2] --> roughly(.2)

	zstar, ustar = suppFcn([-1,  1], w, :Min)
	@fact zstar --> roughly(-.35)
	@fact ustar[1] --> roughly(.1)
	@fact ustar[2] --> roughly(-.25)

	zstar, ustar = suppFcn([-1,  1], w, :Max)
	@fact zstar --> roughly(.25)
	@fact ustar[1] --> roughly(-.05)
	@fact ustar[2] --> roughly(.2)
end

facts("portTest for UM") do
	srand(8675309); data = randn(500, 2)
	w = UMOracle(data, -3.1*ones(2), 2.9 * ones(2), .1, .2)
	portTest(w, -1.7725987075971736, [1, 0])
end


text:
# Automatically generated using Clang.jl

# Skipping MacroDefinition: cpAssertSoft ( __condition__ , ... ) if ( ! ( __condition__ ) ) { cpMessage ( # __condition__ , __FILE__ , __LINE__ , 1 , 0 , __VA_ARGS__ ) ; abort ( ) ; }
# Skipping MacroDefinition: cpAssertWarn ( __condition__ , ... ) if ( ! ( __condition__ ) ) cpMessage ( # __condition__ , __FILE__ , __LINE__ , 0 , 0 , __VA_ARGS__ )
# Skipping MacroDefinition: cpAssertHard ( __condition__ , ... ) if ( ! ( __condition__ ) ) { cpMessage ( # __condition__ , __FILE__ , __LINE__ , 1 , 1 , __VA_ARGS__ ) ; abort ( ) ; }

const CP_USE_DOUBLES = 1

# Skipping MacroDefinition: INFINITY ( __builtin_inf ( ) )
# Skipping MacroDefinition: CP_PI ( ( cpFloat ) 3.14159265358979323846264338327950288 )

const cpTrue = 1
const cpFalse = 0

# Skipping MacroDefinition: CP_NO_GROUP ( ( cpGroup ) 0 )
# Skipping MacroDefinition: CP_ALL_CATEGORIES ( ~ ( cpBitmask ) 0 )
# Skipping MacroDefinition: CP_WILDCARD_COLLISION_TYPE ( ~ ( cpCollisionType ) 0 )

const CP_BUFFER_BYTES = 32 * 1024
const CP_MAX_CONTACTS_PER_ARBITER = 2

# Skipping MacroDefinition: CP_ARBITER_GET_SHAPES ( __arb__ , __a__ , __b__ ) cpShape * __a__ , * __b__ ; cpArbiterGetShapes ( __arb__ , & __a__ , & __b__ ) ;
# Skipping MacroDefinition: CP_ARBITER_GET_BODIES ( __arb__ , __a__ , __b__ ) cpBody * __a__ , * __b__ ; cpArbiterGetBodies ( __arb__ , & __a__ , & __b__ ) ;

const CP_VERSION_MAJOR = 7
const CP_VERSION_MINOR = 0
const CP_VERSION_RELEASE = 3

# Skipping MacroDefinition: CP_CONVEX_HULL ( __count__ , __verts__ , __count_var__ , __verts_var__ ) cpVect * __verts_var__ = ( cpVect * ) alloca ( __count__ * sizeof ( cpVect ) ) ; int __count_var__ = cpConvexHull ( __count__ , __verts__ , __verts_var__ , NULL , 0.0 ) ;

const cpFloat = Cdouble
const cpHashValue = Csize_t
const cpCollisionID = UInt32
const cpBool = Cuchar
const cpDataPointer = Ptr{Cvoid}
const cpCollisionType = Csize_t
const cpGroup = Csize_t
const cpBitmask = UInt32
const cpTimestamp = UInt32

struct cpArbiterFwd end
struct cpArrayFwd end
struct cpBodyFwd end
struct cpCollisionHandlerFwd end
struct cpConstraintClassFwd end
struct cpConstraintFwd end
struct cpContactFwd end
struct cpContactBufferHeaderFwd end
struct cpHashSetFwd end
struct cpShapeClassFwd end
struct cpShapeFwd end
struct cpSpaceFwd end
struct cpSpatialIndexClassFwd end
struct cpSpatialIndexFwd end
struct cpSplittingPlaneFwd end

struct cpVect
    x::cpFloat
    y::cpFloat
end

struct cpTransform
    a::cpFloat
    b::cpFloat
    c::cpFloat
    d::cpFloat
    tx::cpFloat
    ty::cpFloat
end

struct cpMat2x2
    a::cpFloat
    b::cpFloat
    c::cpFloat
    d::cpFloat
end

struct cpArray
    num::Cint
    max::Cint
    arr::Ptr{Ptr{Cvoid}}
end

const cpHashSet = Cvoid
const cpBodyVelocityFunc = Ptr{Cvoid}
const cpBodyPositionFunc = Ptr{Cvoid}
const cpSpatialIndexDestroyImpl = Ptr{Cvoid}
const cpSpatialIndexCountImpl = Ptr{Cvoid}
const cpSpatialIndexEachImpl = Ptr{Cvoid}
const cpSpatialIndexContainsImpl = Ptr{Cvoid}
const cpSpatialIndexInsertImpl = Ptr{Cvoid}
const cpSpatialIndexRemoveImpl = Ptr{Cvoid}
const cpSpatialIndexReindexImpl = Ptr{Cvoid}
const cpSpatialIndexReindexObjectImpl = Ptr{Cvoid}
const cpSpatialIndexReindexQueryImpl = Ptr{Cvoid}
const cpSpatialIndexQueryImpl = Ptr{Cvoid}
const cpSpatialIndexSegmentQueryImpl = Ptr{Cvoid}

struct cpSpatialIndexClass
    destroy::cpSpatialIndexDestroyImpl
    count::cpSpatialIndexCountImpl
    each::cpSpatialIndexEachImpl
    contains::cpSpatialIndexContainsImpl
    insert::cpSpatialIndexInsertImpl
    remove::cpSpatialIndexRemoveImpl
    reindex::cpSpatialIndexReindexImpl
    reindexObject::cpSpatialIndexReindexObjectImpl
    reindexQuery::cpSpatialIndexReindexQueryImpl
    query::cpSpatialIndexQueryImpl
    segmentQuery::cpSpatialIndexSegmentQueryImpl
end

const cpSpatialIndexBBFunc = Ptr{Cvoid}

struct cpSpatialIndex
    klass::Ptr{cpSpatialIndexClassFwd}
    bbfunc::cpSpatialIndexBBFunc
    staticIndex::Ptr{cpSpatialIndexFwd}
    dynamicIndex::Ptr{cpSpatialIndexFwd}
end

const cpContactBufferHeader = Cvoid
const cpCollisionBeginFunc = Ptr{Cvoid}
const cpCollisionPreSolveFunc = Ptr{Cvoid}
const cpCollisionPostSolveFunc = Ptr{Cvoid}
const cpCollisionSeparateFunc = Ptr{Cvoid}

struct cpCollisionHandler
    typeA::cpCollisionType
    typeB::cpCollisionType
    beginFunc::cpCollisionBeginFunc
    preSolveFunc::cpCollisionPreSolveFunc
    postSolveFunc::cpCollisionPostSolveFunc
    separateFunc::cpCollisionSeparateFunc
    userData::cpDataPointer
end

struct cpSpace
    iterations::Cint
    gravity::cpVect
    damping::cpFloat
    idleSpeedThreshold::cpFloat
    sleepTimeThreshold::cpFloat
    collisionSlop::cpFloat
    collisionBias::cpFloat
    collisionPersistence::cpTimestamp
    userData::cpDataPointer
    stamp::cpTimestamp
    curr_dt::cpFloat
    dynamicBodies::Ptr{cpArrayFwd}
    staticBodies::Ptr{cpArrayFwd}
    rousedBodies::Ptr{cpArrayFwd}
    sleepingComponents::Ptr{cpArrayFwd}
    shapeIDCounter::cpHashValue
    staticShapes::Ptr{cpSpatialIndexFwd}
    dynamicShapes::Ptr{cpSpatialIndexFwd}
    constraints::Ptr{cpArrayFwd}
    arbiters::Ptr{cpArrayFwd}
    contactBuffersHead::Ptr{cpContactBufferHeaderFwd}
    cachedArbiters::Ptr{cpHashSetFwd}
    pooledArbiters::Ptr{cpArrayFwd}
    allocatedBuffers::Ptr{cpArrayFwd}
    locked::UInt32
    usesWildcards::cpBool
    collisionHandlers::Ptr{cpHashSetFwd}
    defaultHandler::cpCollisionHandler
    skipPostStep::cpBool
    postStepCallbacks::Ptr{cpArrayFwd}
    staticBody::Ptr{cpBodyFwd}
end

struct cpBodySleeping
    root::Ptr{cpBodyFwd}
    next::Ptr{cpBodyFwd}
    idleTime::cpFloat
end

struct cpBody
    velocity_func::cpBodyVelocityFunc
    position_func::cpBodyPositionFunc
    m::cpFloat
    m_inv::cpFloat
    i::cpFloat
    i_inv::cpFloat
    cog::cpVect
    p::cpVect
    v::cpVect
    f::cpVect
    a::cpFloat
    w::cpFloat
    t::cpFloat
    transform::cpTransform
    userData::cpDataPointer
    v_bias::cpVect
    w_bias::cpFloat
    space::Ptr{cpSpaceFwd}
    shapeList::Ptr{cpShapeFwd}
    arbiterList::Ptr{cpArbiterFwd}
    constraintList::Ptr{cpConstraintFwd}
    sleeping::cpBodySleeping
end

@cenum cpShapeType::UInt32 begin
    CP_CIRCLE_SHAPE = 0
    CP_SEGMENT_SHAPE = 1
    CP_POLY_SHAPE = 2
    CP_NUM_SHAPES = 3
end

const cpShapeCacheDataImpl = Ptr{Cvoid}
const cpShapeDestroyImpl = Ptr{Cvoid}
const cpShapePointQueryImpl = Ptr{Cvoid}
const cpShapeSegmentQueryImpl = Ptr{Cvoid}

struct cpShapeClass
    type::cpShapeType
    cacheData::cpShapeCacheDataImpl
    destroy::cpShapeDestroyImpl
    pointQuery::cpShapePointQueryImpl
    segmentQuery::cpShapeSegmentQueryImpl
end

struct cpShapeMassInfo
    m::cpFloat
    i::cpFloat
    cog::cpVect
    area::cpFloat
end

struct cpBB
    l::cpFloat
    b::cpFloat
    r::cpFloat
    t::cpFloat
end

struct cpShapeFilter
    group::cpGroup
    categories::cpBitmask
    mask::cpBitmask
end

struct cpShape
    klass::Ptr{cpShapeClassFwd}
    space::Ptr{cpSpaceFwd}
    body::Ptr{cpBodyFwd}
    massInfo::cpShapeMassInfo
    bb::cpBB
    sensor::cpBool
    e::cpFloat
    u::cpFloat
    surfaceV::cpVect
    userData::cpDataPointer
    type::cpCollisionType
    filter::cpShapeFilter
    next::Ptr{cpShapeFwd}
    prev::Ptr{cpShapeFwd}
    hashid::cpHashValue
end

struct cpArbiterThread
    next::Ptr{cpArbiterFwd}
    prev::Ptr{cpArbiterFwd}
end

struct cpContact
    r1::cpVect
    r2::cpVect
    nMass::cpFloat
    tMass::cpFloat
    bounce::cpFloat
    jnAcc::cpFloat
    jtAcc::cpFloat
    jBias::cpFloat
    bias::cpFloat
    hash::cpHashValue
end

@cenum cpArbiterState::UInt32 begin
    CP_ARBITER_STATE_FIRST_COLLISION = 0
    CP_ARBITER_STATE_NORMAL = 1
    CP_ARBITER_STATE_IGNORE = 2
    CP_ARBITER_STATE_CACHED = 3
    CP_ARBITER_STATE_INVALIDATED = 4
end

struct cpArbiter
    e::cpFloat
    u::cpFloat
    surface_vr::cpVect
    data::cpDataPointer
    a::Ptr{cpShapeFwd}
    b::Ptr{cpShapeFwd}
    body_a::Ptr{cpBodyFwd}
    body_b::Ptr{cpBodyFwd}
    thread_a::cpArbiterThread
    thread_b::cpArbiterThread
    count::Cint
    contacts::Ptr{cpContactFwd}
    n::cpVect
    handler::Ptr{cpCollisionHandlerFwd}
    handlerA::Ptr{cpCollisionHandlerFwd}
    handlerB::Ptr{cpCollisionHandlerFwd}
    swapped::cpBool
    stamp::cpTimestamp
    state::cpArbiterState
end

const cpConstraintPreStepImpl = Ptr{Cvoid}
const cpConstraintApplyCachedImpulseImpl = Ptr{Cvoid}
const cpConstraintApplyImpulseImpl = Ptr{Cvoid}
const cpConstraintGetImpulseImpl = Ptr{Cvoid}

struct cpConstraintClass
    preStep::cpConstraintPreStepImpl
    applyCachedImpulse::cpConstraintApplyCachedImpulseImpl
    applyImpulse::cpConstraintApplyImpulseImpl
    getImpulse::cpConstraintGetImpulseImpl
end

const cpConstraintPreSolveFunc = Ptr{Cvoid}
const cpConstraintPostSolveFunc = Ptr{Cvoid}

struct cpConstraint
    klass::Ptr{cpConstraintClassFwd}
    space::Ptr{cpSpaceFwd}
    a::Ptr{cpBodyFwd}
    b::Ptr{cpBodyFwd}
    next_a::Ptr{cpConstraintFwd}
    next_b::Ptr{cpConstraintFwd}
    maxForce::cpFloat
    errorBias::cpFloat
    maxBias::cpFloat
    collideBodies::cpBool
    preSolve::cpConstraintPreSolveFunc
    postSolve::cpConstraintPostSolveFunc
    userData::cpDataPointer
end

struct cpCircleShape
    shape::cpShape
    c::cpVect
    tc::cpVect
    r::cpFloat
end

struct cpSegmentShape
    shape::cpShape
    a::cpVect
    b::cpVect
    n::cpVect
    ta::cpVect
    tb::cpVect
    tn::cpVect
    r::cpFloat
    a_tangent::cpVect
    b_tangent::cpVect
end

struct cpSplittingPlane
    v0::cpVect
    n::cpVect
end

struct cpPolyShape
    shape::cpShape
    r::cpFloat
    count::Cint
    planes::Ptr{cpSplittingPlaneFwd}
    _planes::NTuple{12, cpSplittingPlane}
end

struct cpPinJoint
    constraint::cpConstraint
    anchorA::cpVect
    anchorB::cpVect
    dist::cpFloat
    r1::cpVect
    r2::cpVect
    n::cpVect
    nMass::cpFloat
    jnAcc::cpFloat
    bias::cpFloat
end

struct cpSlideJoint
    constraint::cpConstraint
    anchorA::cpVect
    anchorB::cpVect
    min::cpFloat
    max::cpFloat
    r1::cpVect
    r2::cpVect
    n::cpVect
    nMass::cpFloat
    jnAcc::cpFloat
    bias::cpFloat
end

struct cpPivotJoint
    constraint::cpConstraint
    anchorA::cpVect
    anchorB::cpVect
    r1::cpVect
    r2::cpVect
    k::cpMat2x2
    jAcc::cpVect
    bias::cpVect
end

struct cpGrooveJoint
    constraint::cpConstraint
    grv_n::cpVect
    grv_a::cpVect
    grv_b::cpVect
    anchorB::cpVect
    grv_tn::cpVect
    clamp::cpFloat
    r1::cpVect
    r2::cpVect
    k::cpMat2x2
    jAcc::cpVect
    bias::cpVect
end

const cpDampedSpringForceFunc = Ptr{Cvoid}

struct cpDampedSpring
    constraint::cpConstraint
    anchorA::cpVect
    anchorB::cpVect
    restLength::cpFloat
    stiffness::cpFloat
    damping::cpFloat
    springForceFunc::cpDampedSpringForceFunc
    target_vrn::cpFloat
    v_coef::cpFloat
    r1::cpVect
    r2::cpVect
    nMass::cpFloat
    n::cpVect
    jAcc::cpFloat
end

const cpDampedRotarySpringTorqueFunc = Ptr{Cvoid}

struct cpDampedRotarySpring
    constraint::cpConstraint
    restAngle::cpFloat
    stiffness::cpFloat
    damping::cpFloat
    springTorqueFunc::cpDampedRotarySpringTorqueFunc
    target_wrn::cpFloat
    w_coef::cpFloat
    iSum::cpFloat
    jAcc::cpFloat
end

struct cpRotaryLimitJoint
    constraint::cpConstraint
    min::cpFloat
    max::cpFloat
    iSum::cpFloat
    bias::cpFloat
    jAcc::cpFloat
end

struct cpRatchetJoint
    constraint::cpConstraint
    angle::cpFloat
    phase::cpFloat
    ratchet::cpFloat
    iSum::cpFloat
    bias::cpFloat
    jAcc::cpFloat
end

struct cpGearJoint
    constraint::cpConstraint
    phase::cpFloat
    ratio::cpFloat
    ratio_inv::cpFloat
    iSum::cpFloat
    bias::cpFloat
    jAcc::cpFloat
end

const cpSimpleMotorJoint = Cvoid

struct cpContactPointSetPoints
    pointA::cpVect
    pointB::cpVect
    distance::cpFloat
end

struct cpContactPointSet
    count::Cint
    normal::cpVect
    points::cpContactPointSetPoints
end

const cpSpatialIndexIteratorFunc = Ptr{Cvoid}
const cpSpatialIndexQueryFunc = Ptr{Cvoid}
const cpSpatialIndexSegmentQueryFunc = Ptr{Cvoid}
const cpSpaceHash = Cvoid
const cpBBTree = Cvoid
const cpBBTreeVelocityFunc = Ptr{Cvoid}
const cpSweep1D = Cvoid

@cenum cpBodyType::UInt32 begin
    CP_BODY_TYPE_DYNAMIC = 0
    CP_BODY_TYPE_KINEMATIC = 1
    CP_BODY_TYPE_STATIC = 2
end

const cpBodyShapeIteratorFunc = Ptr{Cvoid}
const cpBodyConstraintIteratorFunc = Ptr{Cvoid}
const cpBodyArbiterIteratorFunc = Ptr{Cvoid}

struct cpPointQueryInfo
    shape::Ptr{cpShapeFwd}
    point::cpVect
    distance::cpFloat
    gradient::cpVect
end

struct cpSegmentQueryInfo
    shape::Ptr{cpShapeFwd}
    point::cpVect
    normal::cpVect
    alpha::cpFloat
end

struct cpSimpleMotor
    constraint::cpConstraint
    rate::cpFloat
    iSum::cpFloat
    jAcc::cpFloat
end

const cpPostStepFunc = Ptr{Cvoid}
const cpSpacePointQueryFunc = Ptr{Cvoid}
const cpSpaceSegmentQueryFunc = Ptr{Cvoid}
const cpSpaceBBQueryFunc = Ptr{Cvoid}
const cpSpaceShapeQueryFunc = Ptr{Cvoid}
const cpSpaceBodyIteratorFunc = Ptr{Cvoid}
const cpSpaceShapeIteratorFunc = Ptr{Cvoid}
const cpSpaceConstraintIteratorFunc = Ptr{Cvoid}

struct cpSpaceDebugColor
    r::Cfloat
    g::Cfloat
    b::Cfloat
    a::Cfloat
end

const cpSpaceDebugDrawCircleImpl = Ptr{Cvoid}
const cpSpaceDebugDrawSegmentImpl = Ptr{Cvoid}
const cpSpaceDebugDrawFatSegmentImpl = Ptr{Cvoid}
const cpSpaceDebugDrawPolygonImpl = Ptr{Cvoid}
const cpSpaceDebugDrawDotImpl = Ptr{Cvoid}
const cpSpaceDebugDrawColorForShapeImpl = Ptr{Cvoid}

@cenum cpSpaceDebugDrawFlags::UInt32 begin
    CP_SPACE_DEBUG_DRAW_SHAPES = 1
    CP_SPACE_DEBUG_DRAW_CONSTRAINTS = 2
    CP_SPACE_DEBUG_DRAW_COLLISION_POINTS = 4
end

struct cpSpaceDebugDrawOptions
    drawCircle::cpSpaceDebugDrawCircleImpl
    drawSegment::cpSpaceDebugDrawSegmentImpl
    drawFatSegment::cpSpaceDebugDrawFatSegmentImpl
    drawPolygon::cpSpaceDebugDrawPolygonImpl
    drawDot::cpSpaceDebugDrawDotImpl
    flags::cpSpaceDebugDrawFlags
    shapeOutlineColor::cpSpaceDebugColor
    colorForShape::cpSpaceDebugDrawColorForShapeImpl
    constraintColor::cpSpaceDebugColor
    collisionPointColor::cpSpaceDebugColor
    data::cpDataPointer
end

const cpSpacePointQueryBlock = Cvoid
const cpSpaceSegmentQueryBlock = Cvoid
const cpSpaceBBQueryBlock = Cvoid
const cpSpaceShapeQueryBlock = Cvoid
const CP_POLY_SHAPE_INLINE_ALLOC = 6

struct cpCollisionInfo
    a::Ptr{cpShapeFwd}
    b::Ptr{cpShapeFwd}
    id::cpCollisionID
    n::cpVect
    count::Cint
    arr::Ptr{cpContactFwd}
end

const cpSpaceArbiterApplyImpulseFunc = Ptr{Cvoid}

struct cpPostStepCallback
    func::cpPostStepFunc
    key::Ptr{Cvoid}
    data::Ptr{Cvoid}
end


text:
module ForwardAutodiff

import Base: exp, sin, cos, tan, +, -, *, /, sqrt, convert, promote_rule, zero

export D

struct Infinitesimal{T <: Number} <: Number
    x::T
    dx::T
end

function convert(::Type{Infinitesimal{T}}, x::T) where {T}
   Infinitesimal(x, zero(x))
end
function convert(::Type{Infinitesimal{T}}, x::Infinitesimal{S}) where {S, T}
   Infinitesimal(T(x.x), T(x.dx))
end
# This is needed according to an error before
function convert(::Type{Infinitesimal{T}}, x::T) where {T <: Number}
    Infinitesimal(x, zero(x))
end
function convert(::Type{Infinitesimal{T}}, x::S) where {T, S <: Number}
    x_as_T = convert(T, x)
    Infinitesimal(x_as_T, zero(x_as_T))
end

function zero(x::Infinitesimal{T}) where T
    Infinitesimal(zero(x.x), zero(x.dx))
end

function promote_rule(::Type{Infinitesimal{T}}, ::Type{Infinitesimal{S}}) where {T,S}
    Infinitesimal{promote_type(T,S)}
end
function promote_rule(::Type{Infinitesimal{T}}, ::Type{S}) where {T, S <: Number}
    Infinitesimal{promote_type(T,S)}
end
function promote_rule(::Type{T}, ::Type{Infinitesimal{S}}) where {T <: Number, S}
    Infinitesimal{promote_type(T,S)}
end
function promote_rule(::Type{S}, ::Type{Infinitesimal{T}}) where {S <: AbstractIrrational, T}
    Infinitesimal{promote_type(S,T)}
end

function extract_derivative(xdx::Infinitesimal)
    return xdx.dx
end
function extract_derivative(xs::Array)
    [extract_derivative(x) for x in xs]
end
function extract_derivative(xs::Tuple)
    convert(Tuple, [extract_derivative(x) for x in xs])
end

"""    D([i], f)

Differential operator; with an optional integer, specifies a partial derivative
with respect to argment `i` (from 1).

Returns a function that computes the derivative of `f`.  NOTE: at the moment,
nested derivatives won't work: `D(D(f))` will not compute a second derivative.
"""
function D(f)
    function df(x)
        xdx = Infinitesimal(x, one(x))
        result = f(xdx)
        return extract_derivative(result)
    end
    df
end
function D(i::Integer, f)
    function df(xs...)
        xarr = [(j != i ? x : Infinitesimal(x, one(x))) for (j,x) in enumerate(xs)]
        result = f(xarr...)
        return extract_derivative(result)
    end
    df
end

function +(x::Infinitesimal, y::Infinitesimal)
    Infinitesimal(x.x+y.x, x.dx+y.dx)
end
function -(x::Infinitesimal, y::Infinitesimal)
    Infinitesimal(x.x-y.x, x.dx-y.dx)
end
function -(x::Infinitesimal)
    Infinitesimal(-x.x, -x.dx)
end
function *(x::Infinitesimal, y::Infinitesimal)
    Infinitesimal(x.x*y.x, x.x*y.dx + x.dx*y.x)
end
function /(x::Infinitesimal, y::Infinitesimal)
    Infinitesimal(x.x/y.x, x.dx/y.x - x.x*y.dx/(y.x*y.x))
end

function sqrt(x::Infinitesimal)
    Infinitesimal(sqrt(x.x), x.dx/(2*sqrt(x.x)))
end

function exp(x::Infinitesimal)
    return Infinitesimal(exp(x.x), exp(x.x)*x.dx)
end
function sin(x::Infinitesimal)
    return Infinitesimal(sin(x.x), cos(x.x)*x.dx)
end
function cos(x::Infinitesimal)
    return Infinitesimal(cos(x.x), -sin(x.x)*x.dx)
end
function tan(x::Infinitesimal)
    c = cos(x.x)
    return Infinitesimal(tan(x.x), x.dx/(c*c))
end

end # module


text:
"""
    initialize_matrices(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}})

Initialize the matrices for distances and pathes.

# Arguments
- `graph::Dict{Int64,Array{Tuple{Int64,Int64},1}}`: Graph of the connected nodes with weights
"""
function initialize_matrices(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}})
    # Define intial infinity weight & next matrix for given data type
    vertex_size = find_global_maximum_complex(graph)
    dist = fill(99999, (vertex_size, vertex_size))
    next = zeros(Int64, (vertex_size, vertex_size))

    # Updating the matrix with actual vertex values
    for (key, elements) in graph
        # Setting the diagonal-elements of the weight matrix to 0
        dist[key, key] = 0
        # Setting the diagonal-elements of the next matrix to vertex
        next[key, key] = key
        # Weigthing of the edges and pathes
        for value in elements
            dist[key, value[1]] = value[2]
            next[key, value[1]] = value[1]
        end
    end
    return dist, next, vertex_size
end

"""
    path_reconstruction(next::Array{Int64,2}, u::Int64, v::Int64)

Reconstruction of the actual path between any two endpoint vertices (u & v).

# Arguments
- `graph::Dict{Int64,Array{Tuple{Int64,Int64},1}}`: Graph of the connected nodes with weights
- `u::Int64`: Startpoint of the to investigated path
- `v::Int64`: Endpoint of the to investigated path
"""
function path_reconstruction(next::Array{Int64,2}, u::Int64, v::Int64)
    # Check if path exists
    if next[u, v] == 0
        return zeros(Int64, false)
    end

    # Initialize path
    path = [u]

    while u != v
        u = next[u, v]
        append!(path, u)
    end
    return path
end

"""
    shortest_path_tree(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}}, u=nothing, v=nothing)

The Shortest Path Tree (SPT) algorithm solves the shortest path problem between every
pair of vertices in a given edge-weighted directed Graph based on the Floyd–Warshall
algorithm. Optional, the SPT also provides the total parts between a start- (u) and end-
point (v). For more information see: [https://en.wikipedia.org/wiki/Floyd–Warshall_algorithm](https://en.wikipedia.org/wiki/Floyd–Warshall_algorithm)

# Arguments
- `next::Array{Int64,2}`: Vertex matrix of the connected nodes
- `u::Int64`: Startpoint of the to investigated path; optional
- `v::Int64`: Endpoint of the to investigated path; optional

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> graph = Dict( 1 => [(3, -2)], 2 => [(1, 4),(3, 3)], 3 => [(4, 2)], 4 => [(2, -1)]
julia> ClassicAlgorithmsCollections.shortest_path_tree(graph, 2, 4)
([0 -1 -2 0; 4 0 2 4; 5 1 0 2; 3 -1 1 0], [2, 1, 3, 4])
```
"""
function shortest_path_tree(
    graph::Dict{Int64,Array{Tuple{Int64,Int64},1}},
    u = nothing,
    v = nothing,
)

    dist, next, vertex_size = initialize_matrices(graph)

    # Floyd–Warshall-Algorithm
    for k in 1:vertex_size
        for i in 1:vertex_size
            for j in 1:vertex_size
                # Checking for minimum
                if dist[i, j] > dist[i, k] + dist[k, j]
                    # Add new distance
                    dist[i, j] = dist[i, k] + dist[k, j]
                    # Add new path points
                    next[i, j] = next[i, k]
                end
            end
        end
    end

    # Path-Reconstruction
    if !(isnothing(u) && isnothing(v))
        return dist, path_reconstruction(next, u, v)
    else
        return dist
    end
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using Learn
using Documenter

"""
    doctest_setup()

Importing modules and defining helper functions for doctest
"""
function doctest_setup()
quote
    using Learn
    # Extract the string representation of the method of f having parameters p to
    # make sure the doctest runs correct function
    which_sig(f, p...) = repr(Base.which(f, map(typeof, p))) |> s->split(s," in ") |> first |> strip
end
end

DocMeta.setdocmeta!(Learn, :DocTestSetup, doctest_setup(); recursive=true)

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
# This file is part of TaylorSeries.jl, MIT licensed
#
# Tests for TaylorSeries implementation
using TaylorSeries
using FactCheck
using Compat
import Compat.String

FactCheck.setstyle(:compact)
# FactCheck.onlystats(true)

facts("Tests for Taylor1 expansions") do
    ta(a) = Taylor1([a,one(a)],15)
    t = Taylor1(Int,15)
    tim = im*t
    zt = zero(t)
    ot = 1.0*one(t)
    tol1 = eps(1.0)

    @fact Taylor1([0,1,0,0]) == Taylor1(3)  --> true
    @fact get_coeff(Taylor1(Complex128,3),1) == complex(1.0,0.0) --> true
    @fact eltype(convert(Taylor1{Complex128},ot)) == Complex128  --> true
    @fact eltype(convert(Taylor1{Complex128},1)) == Complex128  --> true
    @fact convert(Taylor1{Complex{Int}},[0,2]) == (2+0im)*t  --> true
    @fact convert(Taylor1{BigFloat},[0.0, 1.0]) == ta(big(0.0))  --> true
    @fact promote(0,Taylor1(1.0,0)) == (zt,ot)  --> true
    @fact eltype(promote(ta(0.0),zeros(Int,2))[2]) == Float64  --> true
    @fact eltype(promote(0,Taylor1(ot))[1]) == Float64  --> true
    @fact eltype(promote(1.0+im, zt)[1]) == Complex{Float64}  --> true
    @fact eltype(TaylorSeries.fixshape(zt,ot)[1]) == Float64  --> true

    @fact length(Taylor1(10)) == 10  --> true
    @fact length(TaylorSeries.fixshape(zt,Taylor1([1.0]))[2]) == 15  --> true
    @fact eltype(TaylorSeries.fixshape(zt,Taylor1([1.0]))[1]) == Float64  --> true
    @fact TaylorSeries.firstnonzero(t) == 1  --> true
    @fact TaylorSeries.firstnonzero(zt) == zt.order+1  --> true

    @fact t == Taylor1(ta(0),15)  --> true
    @fact ot == 1  --> true
    @fact 0.0 == zt  --> true
    @fact get_coeff(tim,1) == complex(0,1)  --> true
    @fact zt+1.0 == ot  --> true
    @fact 1.0-ot == zt  --> true
    @fact t+t == 2t  --> true
    @fact t-t == zt  --> true
    @fact +t == -(-t)  --> true

    tsquare = Taylor1([0,0,1],15)
    @fact t * true == t  --> true
    @fact false * t == zero(t)  --> true
    @fact t^0 == t^0.0 == one(t)  --> true
    @fact t*t == tsquare  --> true
    @fact t*1 == t -->  true
    @fact 0*t == zt -->  true
    @fact (-t)^2 == tsquare  --> true
    @fact t^3 == tsquare*t  --> true
    @fact tsquare/t == t  --> true
    @fact t/(t*3) == (1/3)*ot  --> true
    @fact t/3im == -tim/3  --> true
    @fact 1/(1-t) == Taylor1(ones(t.order+1))  --> true
    @fact Taylor1([0,1,1])/t == t+1  --> true
    @fact (t+im)^2 == tsquare+2im*t-1  --> true
    @fact (t+im)^3 == Taylor1([-1im,-3,3im,1],15)  --> true
    @fact (t+im)^4 == Taylor1([1,-4im,-6,4im,1],15)  --> true
    @fact imag(tsquare+2im*t-1) == 2t  --> true
    @fact (Rational(1,2)*tsquare).coeffs[3] == 1//2  --> true
    @fact t^2/tsquare == ot  --> true
    @fact ((1+t)^(1/3)).coeffs[3]+1/9 <= tol1  --> true
    @fact 1-tsquare == (1+t)-t*(1+t)  --> true
    @fact (1-tsquare)^2 == (1+t)^2.0 * (1-t)^2.0  --> true
    @fact (sqrt(1+t)).coeffs[3] == -1/8  --> true
    @fact ((1-tsquare)^(1//2))^2 == 1-tsquare  --> true
    @fact ((1-t)^(1//4)).coeffs[15] == -4188908511//549755813888  --> true
    @fact abs(((1+t)^3.2).coeffs[14] + 5.4021062656e-5) < tol1  --> true

    @fact isapprox( rem(4.1 + t,4).coeffs[1], (0.1 + t).coeffs[1] )  --> true
    @fact isapprox( mod(4.1 + t,4).coeffs[1], (0.1 + t).coeffs[1] )  --> true
    @fact isapprox( mod2pi(2pi+0.1+t).coeffs[1],(0.1 + t).coeffs[1])  --> true

    @fact abs(ta(1))  --> ta(1)
    @fact abs(ta(-1.0))  --> -ta(-1.0)

    @fact log(exp(tsquare)) == tsquare  --> true
    @fact exp(log(1-tsquare)) == 1-tsquare  --> true
    @fact log((1-t)^2) == 2*log(1-t)  --> true
    @fact real(exp(tim)) == cos(t)  --> true
    @fact imag(exp(tim)) == sin(t)  --> true
    @fact exp(conj(tim)) == cos(t)-im*sin(t) == exp(tim')  --> true
    @fact (exp(t))^(2im) == cos(2t)+im*sin(2t)  --> true
    @fact (exp(t))^Taylor1([-5.2im]) == cos(5.2t)-im*sin(5.2t)  --> true
    @fact get_coeff(convert(Taylor1{Rational{Int}},cos(t)),8) ==
        1//factorial(8)  --> true
    @fact abs((tan(t)).coeffs[8]- 17/315) < tol1  --> true
    @fact abs((tan(t)).coeffs[14]- 21844/6081075) < tol1  --> true
    @fact evaluate(exp(Taylor1([0,1],17)),1.0) == 1.0*e  --> true
    @fact evaluate(exp(Taylor1([0,1],1))) == 1.0  --> true
    @fact evaluate(exp(t),t^2) == exp(t^2)  --> true

    @fact sin(asin(tsquare)) == tsquare --> true
    @fact tan(atan(tsquare)) == tsquare --> true
    @fact atan(tan(tsquare)) == tsquare --> true
    @fact asin(t) + acos(t) == pi/2 --> true
    @fact derivative(acos(t)) == - 1/sqrt(1-t^2) --> true

    @fact - sinh(t) + cosh(t) == exp(-t) --> true
    @fact  sinh(t) + cosh(t) == exp(t) --> true
    @fact evaluate(- sinh(t)^2 + cosh(t)^2 , rand()) == 1 --> true #works in an eps() difference, but ideally it should return a 1 + 𝒪(t¹⁶) Taylor.
    @fact tanh(t + 0im) == -1im * tan(t*1im) --> true
    @fact  evaluate(tanh(t/2),1.5) == evaluate(sinh(t) / (cosh(t) + 1),1.5) --> true
    @fact cosh(t) == real(cos(im*t)) --> true
    @fact sinh(t) == imag(sin(im*t)) --> true

    v = [sin(t), exp(-t)]
    @fact evaluate(v) == [0.0, 1.0]  --> true
    @fact evaluate(v, complex(0.0,0.2)) ==
        [complex(0.0,sinh(0.2)),complex(cos(0.2),sin(-0.2))]  --> true
    @fact derivative(5, exp(ta(1.0))) == exp(1.0)  --> true
    @fact derivative(3, exp(ta(1.0pi))) == exp(1.0pi)  --> true
    @fact isapprox(derivative(10, exp(ta(1.0pi))) , exp(1.0pi) )  --> true
    @fact integrate(derivative(exp(t)),1) == exp(t)  --> true
    @fact integrate(cos(t)) == sin(t)  --> true

    @fact promote(ta(0.0), t) == (ta(0.0),ta(0.0))  --> true

    @fact_throws ArgumentError abs(ta(big(0)))
    @fact_throws ArgumentError 1/t
    @fact_throws ArgumentError zt/zt
    @fact_throws ArgumentError t^1.5
    @fact_throws DomainError t^(-2)
    @fact_throws ArgumentError sqrt(t)
    @fact_throws ArgumentError log(t)
    @fact_throws ArgumentError cos(t)/sin(t)
    @fact_throws AssertionError derivative(30, exp(ta(1.0pi)))

    @fact string(ta(-3)) == " - 3 + 1 t + 𝒪(t¹⁶)"  --> true
    @fact TaylorSeries.pretty_print(ta(3im)) ==
        " ( 3 im )  + ( 1 ) t + 𝒪(t¹⁶)"  --> true
end

facts("Matrix multiplication for Taylor1") do
    order = 30
    n1 = 100
    k1 = 90

    order = max(n1,k1)
    B1 = randn(n1,order)
    Y1 = randn(k1,order)

    A1  = randn(k1,n1)

    for A in (A1,sparse(A1))
        # B and Y contain elements of different orders
        B  = Taylor1{Float64}[Taylor1(collect(B1[i,1:i]),i) for i=1:n1]
        Y  = Taylor1{Float64}[Taylor1(collect(Y1[k,1:k]),k) for k=1:k1]
        Bcopy = deepcopy(B)
        A_mul_B!(Y,A,B)

        # do we get the same result when using the `A*B` form?
        @fact A*B==Y -->true
        # Y should be extended after the multilpication
        @fact reduce(&, [y1.order for y1 in Y] .== Y[1].order) --> true
        # B should be unchanged
        @fact B==Bcopy --> true

        # is the result compatible with the matrix multiplication?  We
        # only check the zeroth order of the Taylor series.
        y1=sum(Y).coeffs[1]
        Y=A*B1[:,1]
        y2=sum(Y)

        # There is a small numerical error when comparing the generic
        # multiplication and the specialized version
        @fact abs(y1-y2) < n1*(eps(y1)+eps(y2)) --> true

        @fact_throws DimensionMismatch A_mul_B!(Y,A[:,1:end-1],B)
        @fact_throws DimensionMismatch A_mul_B!(Y,A[1:end-1,:],B)
        @fact_throws DimensionMismatch A_mul_B!(Y,A,B[1:end-1])
        @fact_throws DimensionMismatch A_mul_B!(Y[1:end-1],A,B)
    end
end

facts("Tests for HomogeneousPolynomial and TaylorN") do

    @fact eltype(set_variables(Int, "x", numvars=2, order=6))  --> TaylorN{Int}
    @fact eltype(set_variables("x", numvars=2, order=6))  --> TaylorN{Float64}
    @fact eltype(set_variables(BigInt, "x y", order=6))  --> TaylorN{BigInt}
    @fact eltype(set_variables("x y", order=6))  --> TaylorN{Float64}

    @fact TaylorSeries.coeff_table[2][1] == [1,0]  --> true
    @fact TaylorSeries.index_table[2][1] == 7  --> true
    @fact TaylorSeries.in_base(get_order(),[2,1]) == 15  --> true
    @fact TaylorSeries.pos_table[4][15] == 2  --> true

    @fact get_order() == 6  --> true
    @fact get_numvars() == 2  --> true

    x, y = set_variables("x y", order=6)
    @fact x.order == 6 --> true
    @fact TaylorSeries.set_variable_names(["x","y"]) == ["x", "y"]  --> true
    @fact TaylorSeries.get_variable_names() == ["x", "y"]  --> true

    set_variables("x", numvars=2, order=17)

    xH = HomogeneousPolynomial([1,0])
    yH = HomogeneousPolynomial([0,1],1)
    xT = TaylorN(xH,17)
    # yT = taylorN_variable(Int64, 2, 17)
    yT = TaylorN(Int64, 2, order=17)
    zeroT = zero( TaylorN([xH],1) )
    uT = one(convert(TaylorN{Float64},yT))
    @fact ones(xH,1) == [HomogeneousPolynomial(1), xH+yH]  --> true
    @fact ones(HomogeneousPolynomial{Complex{Int}},0) ==
        [HomogeneousPolynomial(1+0im)]  --> true
    @fact get_order(zeroT) == 1  --> true
    @fact get_coeff(xT,[1,0]) == 1  --> true
    @fact get_coeff(yH,[1,0]) == 0  --> true
    @fact convert(HomogeneousPolynomial{Int64},[1,1]) == xH+yH  --> true
    @fact convert(HomogeneousPolynomial{Float64},[2,-1]) == 2.0xH-yH  --> true
    @fact convert(TaylorN{Float64}, yH) == 1.0*yT  --> true
    @fact convert(TaylorN{Float64}, [xH,yH]) == xT+1.0*yT  --> true
    @fact convert(TaylorN{Int}, [xH,yH]) == xT+yT  --> true
    @fact promote(xH, [1,1])[2] == xH+yH  --> true
    @fact promote(xH, yT)[1] == xT  --> true
    @fact promote(xT, [xH,yH])[2] == xT+yT  --> true
    @fact typeof(promote(im*xT,[xH,yH])[2]) == TaylorN{Complex{Int64}}  --> true
    # @fact TaylorSeries.fixorder(taylorN_variable(1,1),17) == xT  --> true
    @fact TaylorSeries.fixorder(TaylorN(1, order=1),17) == xT  --> true
    @fact TaylorSeries.iszero(zeroT.coeffs[2])  -->  true

    @fact HomogeneousPolynomial(xH,1) == HomogeneousPolynomial(xH)  --> true
    @fact eltype(xH) == Int  --> true
    @fact length(xH) == 2  --> true
    @fact zero(xH) == 0*xH  --> true
    @fact one(yH) == xH+yH  --> true
    @fact xH * true == xH  --> true
    @fact false * yH == zero(yH)  --> true
    @fact get_order(yH) == 1  --> true
    @fact get_order(xT) == 17  --> true
    @fact xT * true == xT  --> true
    @fact false * yT == zero(yT)  --> true

    @fact xT == TaylorN([xH])  --> true
    @fact one(xT) == TaylorN(1,5)  --> true
    @fact TaylorN(zeroT,5) == 0  --> true
    @fact TaylorN(uT) == convert(TaylorN{Complex},1)  --> true
    @fact get_numvars() == 2  --> true
    @fact length(uT) == get_order()+1  --> true
    @fact eltype(convert(TaylorN{Complex128},1)) == Complex128  --> true

    @fact 1+xT+yT == TaylorN(1,1) + TaylorN([xH,yH],1)  --> true
    @fact xT-yT-1 == TaylorN([-1,xH-yH])  --> true
    @fact xT*yT == TaylorN([HomogeneousPolynomial([0,1,0],2)])  --> true
    @fact (1/(1-xT)).coeffs[4] == HomogeneousPolynomial(1.0,3)  --> true
    @fact xH^20 == HomogeneousPolynomial([0],get_order())  --> true
    @fact (yT/(1-xT)).coeffs[5] == xH^3 * yH  --> true
    @fact mod(1+xT,1) == +xT  --> true
    @fact (rem(1+xT,1)).coeffs[1] == 0  --> true
    @fact abs(1-xT)  --> 1-xT
    @fact abs(-1-xT)  --> 1+xT
    @fact derivative(mod2pi(2pi+yT^3),2) == derivative(yT^3,2)  --> true
    @fact derivative(yT) == zeroT  --> true
    @fact -xT/3im == im*xT/3  --> true
    @fact (xH/3im)' == im*xH/3  --> true

    @fact derivative(2xT*yT^2,1) == 2yT^2  --> true
    @fact xT*xT^3 == xT^4  --> true
    txy = 1.0 + xT*yT - 0.5*xT^2*yT + (1/3)*xT^3*yT + 0.5*xT^2*yT^2
    # @fact (1+taylorN_variable(1,4))^taylorN_variable(2,4) == txy  --> true
    @fact (1+TaylorN(1,order=4))^TaylorN(2,order=4) == txy  --> true
    @fact_throws DomainError yT^(-2)
    @fact_throws DomainError yT^(-2.0)
    @fact (1+xT)^(3//2) == ((1+xT)^0.5)^3  --> true
    @fact real(xH) == xH  --> true
    @fact imag(xH) == zero(xH)  --> true
    @fact conj(im*yH) == (im*yH)'  --> true
    @fact conj(im*yT) == (im*yT)'  --> true
    @fact real( exp(1im * xT)) == cos(xT)  --> true
    @fact get_coeff(convert(TaylorN{Rational{Int}},cos(xT)),[4,0]) ==
        1//factorial(4)  --> true
    cr = convert(TaylorN{Rational{Int}},cos(xT))
    @fact get_coeff(cr,[4,0]) == 1//factorial(4)  --> true
    @fact imag((exp(yT))^(-1im)') == sin(yT)  --> true
    exy = exp( xT+yT )
    @fact evaluate(exy) == 1  --> true
    @fact evaluate(exy,[0.1im,0.01im]) == exp(0.11im)  --> true
    @fact isapprox(evaluate(exy, [1,1]), e^2)  --> true
    txy = tan(xT+yT)
    @fact get_coeff(txy,[8,7]) == 929569/99225  --> true
    ptxy = xT + yT + (1/3)*( xT^3 + yT^3 ) + xT^2*yT + xT*yT^2
    # @fact tan(taylorN_variable(1,4)+taylorN_variable(2,4)) == ptxy  --> true
    @fact tan(TaylorN(1,order=4)+TaylorN(2,order=4)) == ptxy  --> true
    @fact evaluate(xH*yH,[1.0,2.0]) == 2.0  --> true

    g1(xT,yT) = xT^3 + 3yT^2 - 2xT^2 * yT - 7xT + 2
    g2(xT,yT) = yT + xT^2 - xT^4
    f1 = g1(xT,yT)
    f2 = g2(xT,yT)
    @fact gradient(f1) == [ 3*xT^2-4*xT*yT-TaylorN(7,0), 6*yT-2*xT^2 ]  --> true
    @fact ∇(f2) == [2*xT - 4*xT^3, TaylorN(1,0)]  --> true
    @fact jacobian([f1,f2], [2,1]) ==
        jacobian( [g1(xT+2,yT+1), g2(xT+2,yT+1)] )  --> true
    @fact [xT yT]*hessian(f1*f2)*[xT, yT] ==
        [ 2*TaylorN((f1*f2).coeffs[3]) ]  --> true
    @fact hessian(f1^2)/2 == [ [49,0] [0,12] ]  --> true
    @fact hessian(f1-f2-2*f1*f2) == (hessian(f1-f2-2*f1*f2))'  --> true
    @fact hessian(f1-f2,[1,-1]) == hessian(g1(xT+1,yT-1)-g2(xT+1,yT-1))  --> true

    @fact string(-xH) == " - 1 x₁"  --> true
    @fact string(xT^2) == " 1 x₁² + 𝒪(‖x‖¹⁸)"  --> true
    @fact string(1im*yT) == " ( 1 im ) x₂ + 𝒪(‖x‖¹⁸)"  --> true
    @fact string(xT-im*yT) == "  ( 1 ) x₁ - ( 1 im ) x₂ + 𝒪(‖x‖¹⁸)"  --> true

    @fact_throws ArgumentError abs(xT)
    @fact_throws AssertionError 1/x
    @fact_throws AssertionError zero(x)/zero(x)
    @fact_throws ArgumentError sqrt(x)
    @fact_throws AssertionError x^(-2)
    @fact_throws ArgumentError log(x)
    @fact_throws AssertionError cos(x)/sin(y)
end

facts("Testing an identity proved by Euler (8 variables)") do
    make_variable(name, index::Int) = string(name, TaylorSeries.subscriptify(index))

    @compat variable_names = String[make_variable("α", i) for i in 1:4]
    append!(variable_names, [make_variable("β", i) for i in 1:4])

    a1, a2, a3, a4, b1, b2, b3, b4 = set_variables(variable_names, order=4)

    lhs1 = a1^2 + a2^2 + a3^2 + a4^2
    lhs2 = b1^2 + b2^2 + b3^2 + b4^2
    lhs = lhs1 * lhs2

    rhs1 = (a1*b1 - a2*b2 - a3*b3 - a4*b4)^2
    rhs2 = (a1*b2 + a2*b1 + a3*b4 - a4*b3)^2
    rhs3 = (a1*b3 - a2*b4 + a3*b1 + a4*b2)^2
    rhs4 = (a1*b4 + a2*b3 - a3*b2 + a4*b1)^2

    rhs = rhs1 + rhs2 + rhs3 + rhs4

    @fact lhs == rhs  --> true
    v = randn(8)
    @fact evaluate( rhs, v) == evaluate( lhs, v)  --> true
end

facts("High order polynomials test inspired by Fateman (takes a few seconds))") do
    x, y, z, w = set_variables(Int128, "x", numvars=4, order=40)

    function fateman2(degree::Int)
        T = Int128
        oneH = HomogeneousPolynomial(one(T), 0)
        # s = 1 + x + y + z + w
        s = TaylorN(
         [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree)
        s = s^degree
        # s is converted to order 2*ndeg
        s = TaylorN(s, 2*degree)
        return s^2 + s
    end

    function fateman3(degree::Int)
        s = x + y + z + w + 1
        s = s^degree
        s * (s+1)
    end

    f2 = fateman2(20)
    f3 = fateman3(20)
    c = get_coeff(f2,[1,6,7,20])
    @fact c == 128358585324486316800 --> true
    @fact get_coeff(f2,[1,6,7,20]) == c --> true
end

exitstatus()


text:
const SWIPL_LIB = "libswipl"

"""
    Flags

#ifdef PL_KERNEL
#define PL_Q_DEBUG              0x0001  /* = TRUE for backward compatibility */
#endif
#define PL_Q_NORMAL             0x0002  /* normal usage */
#define PL_Q_NODEBUG            0x0004  /* use this one */
#define PL_Q_CATCH_EXCEPTION    0x0008  /* handle exceptions in C */
#define PL_Q_PASS_EXCEPTION     0x0010  /* pass to parent environment */
#define PL_Q_ALLOW_YIELD        0x0020  /* Support I_YIELD */
#define PL_Q_EXT_STATUS         0x0040  /* Return extended status */
#ifdef PL_KERNEL
#define PL_Q_DETERMINISTIC      0x0100  /* call was deterministic */
#endif

"""

const PL_Q_DEBUG = 0x0001
const PL_Q_NORMAL = 0x0002
const PL_Q_NODEBUG = 0x0004
const PL_Q_CATCH_EXCEPTION = 0x0008
const PL_Q_PASS_EXCEPTION = 0x0010
const PL_Q_ALLOW_YIELD = 0x0020
const PL_Q_EXT_STATUS = 0x0040

"""
    Term types

#define PL_VARIABLE      (1)            /* nothing */
#define PL_ATOM          (2)            /* const char * */
#define PL_INTEGER       (3)            /* int */
#define PL_RATIONAL      (4)            /* rational number */
#define PL_FLOAT         (5)            /* double */
#define PL_STRING        (6)            /* const char * */
#define PL_TERM          (7)

#define PL_NIL           (8)            /* The constant [] */
#define PL_BLOB          (9)            /* non-atom blob */
#define PL_LIST_PAIR     (10)           /* [_|_] term */

                                        /* PL_unify_term() */
#define PL_FUNCTOR       (11)           /* functor_t, arg ... */
#define PL_LIST          (12)           /* length, arg ... */
#define PL_CHARS         (13)           /* const char * */
#define PL_POINTER       (14)           /* void * */

"""

const PL_VARIABLE = 1
const PL_ATOM = 2
const PL_INTEGER = 3
const PL_FLOAT = 5
const PL_STRING = 6
const PL_TERM = 7
const PL_NIL = 8
const PL_LIST_PAIR = 10
const PL_FUNCTOR = 11

"""
 Functions for starting and shutting down the SWIPL instance
"""
function swipl_start(exec_path::String)
    argv = [exec_path, "-q"]
    ccall((:PL_initialise, "libswipl"), Cint, (Cint, Ptr{Ptr{UInt8}}), length(argv), argv)
end

swipl_start() = swipl_start("/usr/local/bin/swipl")

function swipl_is_initialised(exec_path::String)
    argv = [exec_path, "-q"]
    ccall((:PL_is_initialised, "libswipl"), Cint, (Cint, Ptr{Ptr{UInt8}}), length(argv), argv)
end

swipl_is_initialised() = swipl_is_initialised()

function swipl_stop(status::Int64)
    ccall((:PL_halt, "libswipl"), Cint, (Cint,), convert(Cint, status))
end

swipl_stop() = swipl_stop(0)

function swipl_cleanup(status::Int64)
    ccall((:PL_cleanup, "libswipl"), Cint, (Cint,), status)
end

swipl_cleanup() = swipl_cleanup(0)

"""
    Atoms and functors
"""

function swipl_new_atom(name::String)
    name_to_use = collect([convert(UInt8, c) for c in name])
    ccall((:PL_new_atom, "libswipl"), Cint, (Cstring,), name) #Cstring(pointer(name)))
end

function swipl_atom_chars(term::Cint)
    r = ccall((:PL_atom_chars, "libswipl"), Cstring, (Cint,), term)
    unsafe_string(r)
end

function swipl_new_functor(atom::Cint, arity::Int64)
    ccall((:PL_new_functor, "libswipl"), Cint, (Cint, Cint), atom, convert(Cint, arity))
end

function swipl_functor_name(functor::Cint)
    ccall((:PL_functor_name, "libswipl"), Cint, (Cint,), functor)
end

function swipl_functor_arity(functor::Cint)
    ccall((:PL_functor_arity, "libswipl"), Cint, (Cint,), functor)
end

"""
    Testing types of terms
"""

function swipl_term_type(term::Cint)
    ccall((:PL_term_type, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_variable(term::Cint)
    ccall((:PL_is_variable, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_ground(term::Cint)
    ccall((:PL_is_ground, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_atom(term::Cint)
    ccall((:PL_is_atom, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_string(term::Cint)
    ccall((:PL_is_string, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_integer(term::Cint)
    ccall((:PL_is_integer, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_float(term::Cint)
    ccall((:PL_is_float, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_callable(term::Cint)
    ccall((:PL_is_callable, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_compound(term::Cint)
    ccall((:PL_is_compound, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_list(term::Cint)
    ccall((:PL_is_list, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_pair(term::Cint)
    ccall((:PL_is_pair, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_atomic(term::Cint)
    ccall((:PL_is_atomic, "libswipl"), Cint, (Cint,), term) > 0
end

function swipl_is_number(term::Cint)
    ccall((:PL_is_number, "libswipl"), Cint, (Cint,), term) > 0
end

"""
    Getting data from terms
"""

function swipl_get_atom(term::Cint)
    atm = Ref{Cint}(0)
    rc = ccall((:PL_get_atom, "libswipl"), Cint, (Cint, Ref{Cint}), term, atm)
    atm[]
end

function swipl_get_integer(term::Cint)
    rint = Ref{Cint}(0)
    ccall((:PL_get_integer, "libswipl"), Cint, (Cint, Ref{Cint}), term, rint)
    rint[]
end

function swipl_get_bool(term::Cint)
    bval = Ref{Cint}()
    ccall((:PL_get_bool, "libswipl"), Cint, (Cint, Ref{Cint}), term, bval)
    bval[]
end

function swipl_get_atom_chrs(term::Cint)
    tmp = Ref{Cstring}(C_NULL)
    rc = ccall((:PL_get_atom_chars, "libswipl"), Cint, (Cint, Ref{Cstring}), term, tmp)
    unsafe_string(tmp[])
end

function swipl_get_float(term::Cint)
    val = Ref{Cfloat}(0)
    ccall((:PL_get_float, "libswipl"), Cint, (Cint, Ref{Cfloat}), term, val)
    val[]
end

function swipl_get_name_arity(term::Cint)
    name = Ref{Cint}(0)
    ar = Ref{Csize_t}(0)
    ccall((:PL_get_name_arity, "libswipl"), Cint, (Cint, Ref{Cint}, Ref{Csize_t}), term, name, ar)
    (name[], ar[])
end

function swipl_get_compound_name_arity(term::Cint)
    name = Ref{Cint}(0)
    ar = Ref{Csize_t}(0)
    ccall((:PL_get_compound_name_arity, SWIPL_LIB), Cint, (Cint, Ref{Cint}, Ref{Csize_t}), term, name, ar)
    (name[], ar[])
end

function swipl_get_arg(index::Csize_t, term::Cint, term_result::Cint)
    ccall((:PL_get_arg, SWIPL_LIB), Cint, (Csize_t, Cint, Cint), index, term, term_result)
end

"""
    Constructing terms
"""
function swipl_new_term_ref()
    ccall((:PL_new_term_ref, "libswipl"), Cint, ())
end

function swipl_new_term_refs(arity::Int64)
    ccall((:PL_new_term_refs, "libswipl"), Cint, (Cint,), arity)
end

function swipl_put_variable(term::Cint)
    ccall((:PL_put_variable, "libswipl"), Cvoid, (Cint,), term)
end

function swipl_put_atom(term::Cint, atom::Cint)
    ccall((:PL_put_atom, "libswipl"), Cvoid, (Cint, Cint), term, atom)
end

function swipl_put_bool(term::Cint, val::Cint)
    ccall((:PL_put_bool, "libswipl"), Cvoid, (Cint, Cint), term, val)
end

function swipl_put_atom_chars(term::Cint, name::String)
    ccall((:PL_put_atom_chars, "libswipl"), Cint, (Cint, Cstring), term, name)
end

function swipl_put_string_chars(term::Cint, name::String)
    ccall((:PL_put_string_chars, "libswipl"), Cint, (Cint, Cstring), term, name)
end

function swipl_put_integer(term::Cint, val::Int64)
    ccall((:PL_put_integer, "libswipl"), Cint, (Cint, Clong), term, val)
end

function swipl_put_float(term::Cint, val::Float64)
    ccall((:PL_put_float, "libswipl"), Cint, (Cint, Cdouble), term, val)
end

function swipl_put_functor(term::Cint, functor::Cint)
    ccall((:PL_put_functor, "libswipl"), Cint, (Cint, Cint), term, functor )
end

function swipl_put_list(term::Cint)
    ccall((:PL_put_list, "libswipl"), Cint, (Cint,), term)
end

function swipl_put_nil(term::Cint)
    ccall((:PL_put_nil, "libswipl"), Cint, (Cint,), term)
end

function swipl_put_term(term1::Cint, term2::Cint)
    ccall((:PL_put_term, "libswipl"), Cint, (Cint, Cint), term1, term2)
end

function swipl_cons_functor(term::Cint, functor::Cint, compound_term::Cint)
    ccall((:PL_cons_functor_v, "libswipl"), Cint, (Cint, Cint, Cint), term, functor, compound_term)
end

function swipl_cons_list(term::Cint, head::Cint, tail::Cint)
    ccall((:PL_cons_list, "libswipl"), Cint, (Cint, Cint, Cint), term, head, tail)
end

function swipl_copy_term_ref(term_from::Cint)
    ccall((:PL_copy_term_ref, SWIPL_LIB), Cint, (Cint,), term_from)
end

"""
    Unification
"""

function swipl_unify(term1::Cint, term2::Cint)
    ccall((:PL_unify, SWIPL_LIB), Cint, (Cint, Cint), term1, term2)
end

function swipl_unify_atom(term1::Cint, term2::Cint)
    ccall((:PL_unify, SWIPL_LIB), Cint, (Cint, Cint), term1, term2)
end

function swipl_unify_bool(term::Cint, val::Bool)
    ccall((:PL_unify_book, SWIPL_LIB), Cint, (Cint, Cint), term, val)
end

function swipl_unify_atom_chars(term::Cint, val::String)
    ccall((:PL_unify_atom_chars, SWIPL_LIB), Cint, (Cint, Cstring), term, val)
end

function swipl_unify_string_chars(term::Cint, val::String)
    ccall((:PL_unify_string_chars, SWIPL_LIB), Cvoid, (Cint, Cstring), term, val)
end

function swipl_unify_integer(term::Cint, val::Int64)
    ccall((:PL_unify_integer, SWIPL_LIB), Cint, (Cint, Cint), term, val)
end

function swipl_unify_float(term::Cint, val::Float64)
    ccall((:PL_unify_float, SWIPL_LIB), Cint, (Cint, Cdouble), term, val)
end

function swipl_unify_functor(term::Cint, functor::Cint)
    ccall((:PL_unify_functor, SWIPL_LIB), Cint, (Cint, Cint), term, functor)
end

function swipl_unify_compound(term::Cint, functor::Cint)
    ccall((:PL_unify_compound, SWIPL_LIB), Cint, (Cint, Cint), term, functor)
end

function swipl_unify_list(term::Cint, head::Cint, tail::Cint)
    ccall((:PL_unify_list, SWIPL_LIB), Cint, (Cint, Cint, Cint), term, head, tail)
end

function swipl_unify_nil(term::Cint)
    ccall((:PL_unify_nil, SWIPL_LIB), Cint, (Cint, ), term)
end

function swipl_unify_arg(index::Int64, term::Cint, arg_term::Cint)
    ccall((:PL_unify_arg, SWIPL_LIB), Cint, (Cint, Cint, Cint), index, term, arg_term)
end

"""
    Calling Prolog
"""

function swipl_predicate(predicate_name::String, arity::Int64)
    ccall((:PL_predicate, SWIPL_LIB), Ptr{Cvoid}, (Cstring, Cint, Cstring), predicate_name, arity, C_NULL)
end

function swipl_open_query(predicate_t::Ptr{Cvoid}, term::Cint, flag::UInt16)
    ccall((:PL_open_query, SWIPL_LIB), Cint, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Cint), C_NULL, flag, predicate_t, term)
end

swipl_open_query(predicate_t::Ptr{Cvoid}, term::Cint) = swipl_open_query(predicate_t, term, PL_Q_NORMAL)

function swipl_next_solution(qid_t::Cint)
    ccall((:PL_next_solution, SWIPL_LIB), Cint, (Cint, ), qid_t)
end

function swipl_cut_query(qid_t::Cint)
    ccall((:PL_cut_query, SWIPL_LIB), Cint, (Cint, ), qid_t)
end

function swipl_close_query(qid_t::Cint)
    ccall((:PL_close_query, SWIPL_LIB), Cint, (Cint,), qid_t)
end


text:
import .DimensionalData
const DD = DimensionalData

axissymbol(::DD.Dim{s}) where {s} = s
axissymbols(o::DD.DimensionalArray) = map(axissymbol, DD.dims(o))
axisvalues(o::DD.DimensionalArray) = DD.val(o)
values(o::DD.DimensionalArray) = parent(o)

function from_namedtuple(::Type{T}, nt::NamedTuple) where {T <: DD.DimensionalArray}
    return DD.DimArray(nt.values, nt.axes)::T
end

roottype(::Type{T}) where {T<:DD.DimensionalArray} = DD.DimensionalArray


text:
function getTableauCrouzeix()
    fac = 0.5/√3
    a = [[ 0.5+fac 0.0    ]
         [-2.0*fac 0.5+fac]]
    b = [0.5,     0.5    ]
    c = [0.5+fac, 0.5-fac]
    o = 2

    TableauDIRK(:crouzeix, o, a, b, c)
end


text:
using Documenter, TopologyPreprocessing

makedocs(sitename="My Documentation")