text:
# A simple function to read user inputs
function input(prompt)
    println(prompt)
    uinput = readline();
    chomp(uinput)
    uinput = parse(Int,uinput)
end

"""
recodeimage(pathtoimage, n_singularvlas)

This function accepts path to a folder containing images (.png & .jpg) and
    number of singular values you want as arguments and generates a CSV File
    in the same folder with each row corresoponding to each images in the folder
    and each column corresoponding to it's respective singular values.
    USAGE:
    ```julia
    recodeimage(path,3)
    ```
    In the above example, we are asking to return a CSV File that has the first 3 singular
    values of the images in the folder.
"""
function recodeimage(pathtoimage)
    # The following code block looks for .jpg/.JPG/.png/.PNG files and create a list of them
    imagelist = glob("*.jpg", pathtoimage)
    append!(imagelist,  glob("*.JPG", pathtoimage))
    append!(imagelist,  glob("*.png", pathtoimage))
    append!(imagelist,  glob("*.PNG", pathtoimage))

    # Algorithm for scree plot ===================================================================================
    # creating a temporary array of arrays to hold the numerical values of the grayscaled images
    # temp[1] will have the grayscaled information of image 1, and so on...
    temp = Array{Array{Float64,2},1}(undef,length(imagelist))

    # @showprogress is a macro to print the progress of this loop when this function is run
        @showprogress "Loading Images " for i in 1:length(imagelist)
        temp[i] = Float64.(Gray.(imresize(load(imagelist[i]),128,128)))
        # TODO: delete the remaining commented code in this loop.
        # @pipe is a macro for chaining multiple tasks
        # the next two blocks of code takes the image, converts it into grayscale and compute
        # the singular values, then only the first n_singular values are stored in the recodedArray
        #img_singluar = @pipe X |> Float64.(Gray.(_)) |> svdvals(_)[1:n_singularvlas]'
        #recodedArray = vcat(recodedArray, img_singluar)
    end

    # stacking individual images to create a giant image matrix
    stackedX = vcat(temp...)
    println("Running Singular Value Decomposition...")
    U,S,V= svd(stackedX)
    S = S ./norm(S)
    screeplot(S)
    # ============================================================================================================

    # After examining the scree plot, the user decides the no. of singular values
    n_singularvlas = input("No. of Features (due to bug in the code that reads user inputs, you might have to enter the no twice, if the program didn't run first time)")

    # Initializing an empty array to store the n_singular values of the images
    recodedArray = Array{Float64}(undef, 0, n_singularvlas)
    # computing singular values using only the first n_singulars
    for imagearrays in temp
        img_singular = @pipe imagearrays |> svdvals(_)[1:n_singularvlas]'
        recodedArray = vcat(recodedArray, img_singular)
    end

    # Normalizing the singular values
    recodedArray =  eachcol(recodedArray) ./ norm.(eachcol(recodedArray))
    # writing the array as a .csv file
    filename = joinpath(pathtoimage, "image_recoded.csv")
    CSV.write(filename,  DataFrame(recodedArray), writeheader=true)
end

"""
recodeaudio(filepath,ncep)

This function reads audio files (.wav files) in the file path and generates an excel file (.csv) with
    each row corresoponding to a particular audio file and each column corresoponding to the mean
    cepstral coefficient. Only the first ncep coefficients are computed.

    USAGE:
    ```julia
    recodeaudio(path,3)
    ```
"""
function recodeaudio(filepath)

    # creating a list of .wav files in the folder and an empty array to store the mean cepstral coefficients
    audiolist = glob("*.wav",filepath)
    recodedArray = Array{Float64}(undef, 0, 26)

    @showprogress for files in audiolist
        x,fs = wavread(files)

        # mfcc is a function from the MFCC package and computes the Mel Frequency Cepstral Coefficients
        # mfcc returns 3 values: [1] a matrix of numcep columns with for each speech frame a row of MFCC coefficient,
        # [2] the power spectrum, and [3] a dictionary containing the information about the parameters
        # mfcc(x,fs, numcep = ncep)[1]...here the [1] specifies that we accesing only the cep coefficients matrix.
        # From the above step, the mean cep is calculated, then the first 13 deltas are extracted
        # from the signal and mean_cep and deltas are horizontally concatinated to form the feature vector
        deltavals = deltas(x)[1:13]'
        cep = collect(mfcc(x,fs)[1])
        features = [mean(cep,dims=1) deltavals] # mean(_,dims=1) returns the column means
        recodedArray = vcat(recodedArray, features)
    end
    replace!(recodedArray, Inf =>0)
    replace!(recodedArray, -Inf =>0)
    replace!(recodedArray, NaN =>0)
    U,S,V = svd(recodedArray)
    S = S ./norm(S)
    screeplot(S)
    # ============================================================================================================

    # After examining the scree plot, the user decides the no. of singular values
    n_singularvlas = input("No. of Features (due to bug in the code that reads user inputs, you might have to enter the no twice, if the program didn't run first time)")
    # Saving the recodedArray as a .csv file
    featurematrix = U[:,1:n_singularvlas]
    filename = joinpath(filepath, "audio_recoded.csv")

    # Normalizing values before writing onto excel sheet.
    #recodedArray = eachcol(recodedArray) ./ norm(eachcol(recodedArray))
    CSV.write(filename,  DataFrame(featurematrix), writeheader=true)
end

"""
recodetext(pathtotxt, n_singularvals)

This function reads PDF files in the pathtotxt and generates an excel file (.CSV) which contains the
    document by term matrix reconstructed using only the first n_singularvals. If the passed n_singularvals
    is a larger number than the singular values present, all singular values are used.

    USAGE:
    ```julia
    recodetext(pathtotxt,3)
    ```
"""
function recodetext(pathtotxt)

    # PDF parsing utilities for this function are provided by the Taro julia package
    # and Taro requires to initialize a java virtual machine. Taro.init() achieves that
    # Creating a list of PDF files in the path.
    files = glob("*.pdf", pathtotxt)

    # Initializing an empty array to store the text from each document
    docs = Array{StringDocument{String},1}[]

    # This block of code extracts the text from the PDFs and creates a corpus
    getTitle(t) = TextAnalysis.sentence_tokenize(Languages.English(), t)[1]
   @showprogress for i in 1:length(files)
        meta,txt = Taro.extract(files[i]);
        txt = replace(txt, '\n' => ' ')
        title = getTitle(txt)
        dm = TextAnalysis.DocumentMetadata(Languages.English(), title, "", meta["Creation-Date"] )
        doc = StringDocument(txt, dm)
        docs = vcat(docs, doc)
    end
    crps = Corpus(docs)

    # cleaning the corpus. getting rid of non-words, punctuations etc...
    prepare!(crps, strip_non_letters | strip_punctuation | strip_case | strip_stopwords)
    TextAnalysis.stem!(crps)
    update_lexicon!(crps)
    update_inverse_index!(crps)

    # The corpus is used to generates a document term matrix,
    # and we perform singluar value decomposition on this matrix.
    U,S,V = @pipe crps |> DocumentTermMatrix(_) |> dtm(_, :dense) |> svd(_)

    S = S ./norm(S)
    p1 = screeplot(S)
    display(p1)
    n_singularvlas = input("No. of Features (due to bug in the code that reads user inputs, you might have to enter the no twice, if the program didn't run first time)")
    # Reconstructing the document term matrix with just the first n_singularvals.
   # DTM = U[:,1:n_singularvals]*diagm(S[1:n_singularvals])*Vt[:,1:n_singularvals]'
    # saving the reconstructed document term matrix as a .csv file
    filename = joinpath(pathtotxt, "reconstructed_document_term_matrix.csv")
    CSV.write(filename,  DataFrame(U[:,1:n_singularvlas]), writeheader=true)
end

# ===================================
function recodeaudio2(filepath)

    # creating a list of .wav files in the folder and an empty array to store the mean cepstral coefficients
    audiolist = glob("*.wav",filepath)
    recodedArray = Array{Float64}(undef, 0, 26)

    @showprogress for files in audiolist
        x,fs = wavread(files)

        # mfcc is a function from the MFCC package and computes the Mel Frequency Cepstral Coefficients
        # mfcc returns 3 values: [1] a matrix of numcep columns with for each speech frame a row of MFCC coefficient,
        # [2] the power spectrum, and [3] a dictionary containing the information about the parameters
        # mfcc(x,fs, numcep = ncep)[1]...here the [1] specifies that we accesing only the cep coefficients matrix.
        # From the above step, the mean cep is calculated, then the first 13 deltas are extracted
        # from the signal and mean_cep and deltas are horizontally concatinated to form the feature vector
        deltavals = deltas(x)[1:13]'
        cep = collect(mfcc(x,fs)[1])
        features = [mean(cep[1],dims=1) deltavals] # mean(_,dims=1) returns the column means
        recodedArray = vcat(recodedArray, features)
    end
    replace!(recodedArray, Inf =>0)
    replace!(recodedArray, -Inf =>0)
    replace!(recodedArray, NaN =>0)
    U,S,V = svd(recodedArray)
    S = S ./norm(S)
    screeplot(S)
    # ============================================================================================================

    # After examining the scree plot, the user decides the no. of singular values
    n_singularvlas = input("No. of Features (due to bug in the code that reads user inputs, you might have to enter the no twice, if the program didn't run first time)")
    # Saving the recodedArray as a .csv file
    featurematrix = U[:,1:n_singularvlas]
    filename = joinpath(filepath, "audio_recoded.csv")

    # Normalizing values before writing onto excel sheet.
    #recodedArray = eachcol(recodedArray) ./ norm(eachcol(recodedArray))
    CSV.write(filename,  DataFrame(featurematrix), writeheader=true)
end


text:
using DataFrames, Missings, CSV

# constructing a DataFrame:
df = DataFrame()
df[:Col1] = 1:4
df[:Col2] = [exp(1), pi, sqrt(2), 42]
df[:Col3] = [true, false, true, false]
show(df)
# 4x3 DataFrame
# | Row | Col1 | Col2    | Col3  |
# |-----|------|---------|-------|
# | 1   | 1    | 2.71828 | true  |
# | 2   | 2    | 3.14159 | false |
# | 3   | 3    | 1.41421 | true  |
# | 4   | 4    | 42.0    | false |

df = DataFrame(Col1 = 1:4, Col2 = [exp(1), pi, sqrt(2), 42], Col3 = [true, false, true, false])

# columns:
show(df[2])
#> [2.718281828459045,3.141592653589793,1.4142135623730951,42.0]
show(df[:Col2]) #> same output

# rows:
show(df[1, :])  #> 1st row
# 1x3 DataFrame
# | Row | Col1 | Col2    | Col3 |
# |-----|------|---------|------|
# | 1   | 1    | 2.71828 | true |
show(df[2:3, :]) #> 2nd and 3rd row:
# 2x3 DataFrame
# | Row | Col1 | Col2    | Col3  |
# |-----|------|---------|-------|
# | 1   | 2    | 3.14159 | false |
# | 2   | 3    | 1.41421 | true  |
show(df[2:3, :Col2]) #> 2nd and 3rd row, only 2nd column:
# [3.141592653589793,1.4142135623730951]
df[2:3, [:Col2, :Col3]] #> 2nd and 3rd row, 2nd and 3rd column:
# 2x2 DataFrame
# | Row | Col2    | Col3  |
# |-----|---------|-------|
# | 1   | 3.14159 | false |
# | 2   | 1.41421 | true  |
head(df)

df0 = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b",  "c"], 10))
head(df0)

tail(df)

names(df)
# 3-element Array{Symbol,1}:
#  :Col1
#  :Col2
#  :Col3
eltypes(df)
# 3-element Array{Type{T<:Top},1}:
 # Int64
 # Float64
 # Bool

describe(df)
# 3×8 DataFrame
# │ Row │ variable │ mean    │ min     │ median  │ max  │ nunique │ nmissing │ eltype  │
# ├─────┼──────────┼─────────┼─────────┼─────────┼──────┼─────────┼──────────┼─────────┤
# │ 1   │ Col1     │ 2.5     │ 1       │ 2.5     │ 4    │         │          │ Int64   │
# │ 2   │ Col2     │ 12.3185 │ 1.41421 │ 2.92994 │ 42.0 │         │          │ Float64 │
# │ 3   │ Col3     │ 0.5     │ false   │ 0.5     │ true │         │          │ Bool    │

fname = "winequality.csv"
data = CSV.read(fname, delim = ';')
# 1599x12 DataFrame
# | Row  | fixed_acidity | volatile_acidity | citric_acid | residual_sugar |
# |------|---------------|------------------|-------------|----------------|
# | 1    | 7.4           | 0.7              | 0.0         | 1.9            |
# | 2    | 7.8           | 0.88             | 0.0         | 2.6            |
# | 3    | 7.8           | 0.76             | 0.04        | 2.3            |
# ⋮
# | 1596 | 5.9           | 0.55             | 0.1         | 2.2            |
# | 1597 | 6.3           | 0.51             | 0.13        | 2.3            |
# | 1598 | 5.9           | 0.645            | 0.12        | 2.0            |
# | 1599 | 6.0           | 0.31             | 0.47        | 3.6            |

# | Row  | chlorides | free_sulfur_dioxide | total_sulfur_dioxide | density |
# |------|-----------|---------------------|----------------------|---------|
# | 1    | 0.076     | 11.0                | 34.0                 | 0.9978  |
# | 2    | 0.098     | 25.0                | 67.0                 | 0.9968  |
# | 3    | 0.092     | 15.0                | 54.0                 | 0.997   |
# ⋮
# | 1596 | 0.062     | 39.0                | 51.0                 | 0.99512 |
# | 1597 | 0.076     | 29.0                | 40.0                 | 0.99574 |
# | 1598 | 0.075     | 32.0                | 44.0                 | 0.99547 |
# | 1599 | 0.067     | 18.0                | 42.0                 | 0.99549 |

# | Row  | pH   | sulphates | alcohol | quality |
# |------|------|-----------|---------|---------|
# | 1    | 3.51 | 0.56      | 9.4     | 5       |
# | 2    | 3.2  | 0.68      | 9.8     | 5       |
# | 3    | 3.26 | 0.65      | 9.8     | 5       |
# ⋮
# | 1596 | 3.52 | 0.76      | 11.2    | 6       |
# | 1597 | 3.42 | 0.75      | 11.0    | 6       |
# | 1598 | 3.57 | 0.71      | 10.2    | 5       |
# | 1599 | 3.39 | 0.66      | 11.0    | 6       |
typeof(data) # DataFrame
size(data) # (1599,12)
data[:density]
# 1599-element DataArray{Float64,1}:
#  0.9978
#  0.9968
#  0.997
#  0.998
#  ⋮
#  0.99512
#  0.99574
#  0.99547
#  0.99549
CSV.write("dataframe1.csv", df,  delim = ';')

# queries:
# the quality for all wines:
data[:quality]
# give the wines with alcohol % = 9.5
show(data[ data[:alcohol].== 9.5, :])
# count the number of wines grouped by quality
show(by(data, :quality, data -> size(data, 1)))
# 6x2 DataFrame
# | Row | quality | x1  |
# |-----|---------|-----|
# | 1   | 3       | 10  |
# | 2   | 4       | 53  |
# | 3   | 5       | 681 |
# | 4   | 6       | 638 |
# | 5   | 7       | 199 |
# | 6   | 8       | 18  |
_, count = hist(data[:quality])
#>  count 6-element Array{Int64,1}:  10 53 681 638 199 18
class = sort(unique(data[:quality]))
#> 6-element DataArray{Int64,1}: 3 4 5 6 7 8
df_quality = DataFrame(qual=class, no=count)
# 6x2 DataFrame
# | Row | qual | no  |
# |-----|------|-----|
# | 1   | 3    | 10  |
# | 2   | 4    | 53  |
# | 3   | 5    | 681 |
# | 4   | 6    | 638 |
# | 5   | 7    | 199 |
# | 6   | 8    | 18  |

ext: jl


text:
using Documenter
using SPHtoGrid

makedocs(
    sitename = "SPHtoGrid.jl",
    format = Documenter.HTML(),
    modules = [SPHtoGrid],
    pages = [
            "Table of Contents" => "index.md",
            "Install" => "install.md",
            "Mapping SPH Data" => "mapping.md",
            "Effect maps" => "effects.md",
            "Rotating Images" => "rotating.md",
            "Saving/Loading Images" => "io.md",
            "External Programs" => "external.md",
            "API reference" => "api.md"
            ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/LudwigBoess/SPHtoGrid.jl.git"
)


text:
# Test rotation to great circle path
using Test
using SAC

@testset "Great circle path rotation" begin
    let
        # The azimuth and distance of the event
        local az = rand(0:359) # degrees, approximate for now
        local d = 0.1 # degrees

        # Pair of SAC traces representing a spike arriving on the radial component
        local npts = 10
        local delta = 1
        local s = SACtr.([delta,delta], [npts,npts])
        s[:stlo], s[:stla] = 0, 0
        s[:evla] = -d*cosd(az)
        s[:evlo] = -d*sind(az)

        # Make azimuth more accurate as using great circle calculation
        az = mod(s[1][:baz] + 180, 360)
        s[:t] = zeros.(s[:npts])
        # Randomly swap the order of the components, which are in random orientations
        local ne = rand(Bool)
        local az1 = rand(0:359)
        s[:cmpaz] = ne ? [az1, az1+90] : [az1+90, az1]
        s[:kcmpnm] = ne ? ["1", "2"] : ["2", "1"]
        s[:cmpinc] = 90
        local imax = npts÷2
        s[1].t[imax], s[2].t[imax] = ne ? (cosd(az-az1), sind(az-az1)) :
                                          (sind(az-az1), cosd(az-az1))

        # Create all possible combinations of rotations
        local (r, t) = rotate_to_gcp(s[1], s[2])
        local (r′, t′) = rotate_to_gcp(s[2], s[1])
        local rt = rotate_to_gcp(s)
        local rt′ = rotate_to_gcp(s[[2,1]])
        local (R, T) = rotate_to_gcp!(deepcopy(s[1]), deepcopy(s[2]))
        local (R′, T′) = rotate_to_gcp!(deepcopy(s[2]), deepcopy(s[1]))
        local RT = rotate_to_gcp!(deepcopy(s))
        local RT′ = rotate_to_gcp!(deepcopy(s[[2,1]]))

        local list_of_radials = [r, r′, rt[1], rt′[1], R, R′, RT[1], RT′[1]]
        local list_of_transverses = [t, t′, rt[2], rt′[2], T, T′, RT[2], RT′[2]]

        # Test rotation has been done correctly
        local atol = sqrt(eps(SAC.SACFloat))
        for (r1,t1) in zip(list_of_radials, list_of_transverses)
            @test r1.t[imax] ≈ 1
            @test all(isapprox.([r1.t[1:imax-1]; r1.t[imax+1:end]], 0, atol=atol))
            @test all(isapprox.(t1.t, 0, atol=atol))
        end

        # Test order of arguments doesn't stop returns being in order r, t
        for r1 in list_of_radials, r2 in list_of_radials
            @test r1 == r2
        end
        for t1 in list_of_transverses, t2 in list_of_transverses
            @test t1 == t2
        end
    end
end


text:
names = [
    "test_sortedindices.jl",
    "test_polynomials.jl",

    "test_bases.jl",
    "test_states.jl",

    "test_operators.jl",
    "test_operators_dense.jl",
    "test_sparsematrix.jl",
    "test_operators_sparse.jl",
    "test_operators_lazytensor.jl",
    "test_operators_lazysum.jl",
    "test_operators_lazyproduct.jl",

    "test_fock.jl",
    "test_spin.jl",
    "test_particle.jl",
    "test_manybody.jl",
    "test_nlevel.jl",
    "test_subspace.jl",
    "test_state_definitions.jl",

    "test_phasespace.jl",
    "test_transformations.jl",

    "test_metrics.jl",
    "test_embed.jl",
    "test_spectralanalysis.jl",

    "test_timeevolution_schroedinger.jl",
    "test_timeevolution_master.jl",
    "test_timeevolution_mcwf.jl",

    "test_timeevolution_twolevel.jl",
    "test_timeevolution_pumpedcavity.jl",

    "test_superoperators.jl",
    "test_steadystate.jl",
    "test_timecorrelations.jl",

    "test_semiclassical.jl",

    "test_stochastic_definitions.jl",
    "test_stochastic_schroedinger.jl",
    "test_stochastic_master.jl",
    "test_stochastic_semiclassical.jl",

    "test_printing.jl"
]

detected_tests = filter(
    name->startswith(name, "test_") && endswith(name, ".jl"),
    readdir("."))

unused_tests = setdiff(detected_tests, names)
if length(unused_tests) != 0
    error("The following tests are not used:\n", join(unused_tests, "\n"))
end

unavailable_tests = setdiff(names, detected_tests)
if length(unavailable_tests) != 0
    error("The following tests could not be found:\n", join(unavailable_tests, "\n"))
end

for name=names
    if startswith(name, "test_") && endswith(name, ".jl")
        include(name)
    end
end


text:
using GlobalSensitivity, Test

function ishi_batch(X)
    A= 7
    B= 0.1
    @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
end
function ishi(X)
    A= 7
    B= 0.1
    sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
end

function linear_batch(X)
    A= 7
    B= 0.1
    @. A*X[1,:]+B*X[2,:]
end
function linear(X)
    A= 7
    B= 0.1
    A*X[1]+B*X[2]
end

lb = -ones(4)*π
ub = ones(4)*π

res1 = gsa(ishi,EASI(),[[lb[i],ub[i]] for i in 1:4],N=15000)
res2 = gsa(ishi_batch,EASI(),[[lb[i],ub[i]] for i in 1:4],N=15000,batch=true)

res1efast = gsa(ishi,eFAST(),[[lb[i],ub[i]] for i in 1:4],n=15000)
res2efast = gsa(ishi_batch,eFAST(),[[lb[i],ub[i]] for i in 1:4],n=15000,batch=true)

@test res1.S1 ≈ res1efast.S1[1,:] atol = 3e-2
@test res2.S1 ≈ res2efast.S1[1,:] atol = 3e-2

res1 = gsa(linear,EASI(),[[lb[i],ub[i]] for i in 1:4], N=15000)
res2 = gsa(linear_batch,EASI(),[[lb[i],ub[i]] for i in 1:4],batch=true, N=15000)

res1efast = gsa(linear,eFAST(),[[lb[i],ub[i]] for i in 1:4])
res2efast = gsa(linear_batch,eFAST(),[[lb[i],ub[i]] for i in 1:4],batch=true)

@test res1.S1 ≈ res1efast.S1[1,:] atol = 3e-2
@test res2.S1 ≈ res2efast.S1[1,:] atol = 3e-2

# function ishi_linear(X)
#     A= 7
#     B= 0.1
#     [sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1]),A*X[1]+B*X[2]]
# end

# function ishi_linear_batch(X)
#     A= 7
#     B= 0.1
#     X1 = @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
#     X2 = @. A*X[1,:]+B*X[2,:]
#     vcat(X1',X2')
# end

# res1 = gsa(ishi_linear,EASI(),[[lb[i],ub[i]] for i in 1:4],N=15000)
# res2 = gsa(ishi_linear_batch,EASI(),[[lb[i],ub[i]] for i in 1:4],N=15000,batch=true)

# res1 = gsa(ishi_linear,eFAST(),[[lb[i],ub[i]] for i in 1:4],n=15000)
# res2 = gsa(ishi_linear_batch,eFAST(),[[lb[i],ub[i]] for i in 1:4],n=15000,batch=true)


text:
"""
"""
struct SingleFieldFEFunction{T<:CellField} <: CellField
  cell_field::T
  cell_vals
  free_values
  dirichlet_values
  fe_space
  @doc """
  """
  function SingleFieldFEFunction(
    cell_field::CellField,
    cell_vals::AbstractArray{<:AbstractArray},
    free_values::AbstractVector,
    dirichlet_values::AbstractVector,
    fe_space::SingleFieldFESpace)

    T = typeof(cell_field)
    new{T}(cell_field,cell_vals,free_values,dirichlet_values,fe_space)
  end
end

FEFunctionStyle(::Type{<:SingleFieldFEFunction}) = Val{true}()

get_array(f::SingleFieldFEFunction) = get_array(f.cell_field)

get_free_values(f::SingleFieldFEFunction) = f.free_values

get_dirichlet_values(f::SingleFieldFEFunction) = f.dirichlet_values

get_fe_space(f::SingleFieldFEFunction) = f.fe_space

get_cell_map(f::SingleFieldFEFunction) = get_cell_map(f.cell_field)

get_cell_values(f::SingleFieldFEFunction) = f.cell_vals

RefStyle(::Type{SingleFieldFEFunction{T}}) where T = RefStyle(T)

CellData.MetaSizeStyle(::Type{SingleFieldFEFunction{T}}) where T = MetaSizeStyle(T)

CellData.get_cell_axes(a::SingleFieldFEFunction) = get_cell_axes(a.cell_field)

CellData.get_memo(a::SingleFieldFEFunction) = get_memo(a.cell_field)


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

using LaTeXStrings
using Plots
pyplot()

# pick a pixel (xprime, yprime), and trace a full ray through the disk

function plot_ray(xprime, yprime)

    v0 = 1.0 # km/s

    # Get a decent estimate of DeltaVmax
    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # DeltaVmax = 0.5 # km/s
    rmax = 600 * AU
    nzp = 512
    zps = linspace(-600 * AU, 600 * AU, nzp)

    rcyls = Array{Float64}(nzp)
    zs = Array{Float64}(nzp)
    Ss = Array{Float64}(nzp)
    Deltavs = Array{Float64}(nzp)
    tildeDeltavs = Array{Float64}(nzp)
    DeltaVs = Array{Float64}(nzp)
    rhos = Array{Float64}(nzp)
    Us = Array{Float64}(nzp)
    tildeUs = Array{Float64}(nzp)
    Deltas = Array{Float64}(nzp)
    alphas = Array{Float64}(nzp)
    tildealphas = Array{Float64}(nzp)
    tildealphasG = Array{Float64}(nzp)

    z0s = DiskTracer.geometry.get_zps(xprime, yprime, pars, v0)
    println("z0s ", z0s)

    # Get the minimum rcyl
    rcyl_min = xprime

    # Occurs at what value of zprime?
    zprime_rcyl_min = cosd(pars.incl)/sind(pars.incl) * yprime

    # Actually, we want the zprime that corresponds to the midplane crossing
    zprime_midplane = -tand(pars.incl) * yprime
    # Calculate the rcyl_mid here

    # Get amplitude at midplane crossing
    # Corresponds to r_cyl_mid, 0.0
    r_cyl_mid = sqrt(xprime^2 + (cosd(pars.incl) * yprime + sind(pars.incl) * zprime_midplane)^2)

    println("rcyl min ", rcyl_min/AU, " zprime_rcyl_min ", zprime_rcyl_min/AU, " zprime_midplane ", zprime_midplane/AU, " rcyl_mid ", r_cyl_mid/AU)

    # Get scale height at rcyl_min (will be an underestimate)
    H_rcyl_min = DiskTracer.model.Hp(rcyl_min, pars)
    H_rcyl_mid = DiskTracer.model.Hp(r_cyl_mid, pars)

    # The amplitude of tilde Upsilon is just the value at the midplane
    Upsilon_amp = DiskTracer.model.Upsilon_nu(r_cyl_mid, 0.0, mol.nu_0, pars, mol)

    # The sigma of the Gaussian for tilde Upsilon is this lengthened by the inclination
    Upsilon_sigma = H_rcyl_mid / cosd(pars.incl)

    sigma_Upsilon = H_rcyl_mid / cosd(pars.incl)
    a_Upsilon = DiskTracer.model.Upsilon_nu(r_cyl_mid, 0.0, mol.nu_0, pars, mol) * sqrt(2pi) * sigma_Upsilon
    mu_Upsilon = zprime_midplane

    DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)

    mu_upsilon1, mu_upsilon2 = z0s
    sigma_upsilon = DeltaV * sqrt(2)/3 * (xprime * sqrt(G * pars.M_star * M_sun) * sind(pars.incl))^(4./3) / (mu_upsilon1 * (v0 * 1.e5)^(7/3))
    a_upsilon = sqrt(2pi) * sigma_upsilon

    sigma_c2 = 1/(sigma_Upsilon^(-2) + sigma_upsilon^(-2))
    mu_c1 = sigma_c2 * (mu_Upsilon/sigma_Upsilon^2 + mu_upsilon1/sigma_upsilon^2)
    mu_c2 = sigma_c2 * (mu_Upsilon/sigma_Upsilon^2 + mu_upsilon2/sigma_upsilon^2)

    a_alpha1 = a_Upsilon * a_upsilon / (sqrt(2pi) * sqrt(sigma_Upsilon^2 + sigma_upsilon^2)) * exp(- (mu_Upsilon - mu_upsilon1)^2/(2 * (sigma_Upsilon^2 + sigma_upsilon^2)))
    a_alpha2 = a_Upsilon * a_upsilon / (sqrt(2pi) * sqrt(sigma_Upsilon^2 + sigma_upsilon^2)) * exp(- (mu_Upsilon - mu_upsilon2)^2/(2 * (sigma_Upsilon^2 + sigma_upsilon^2)))

    for i=1:nzp
        zprime = zps[i]
        rcyl = DiskTracer.geometry.get_r_cyl(xprime, yprime, zprime, pars.incl)
        z = DiskTracer.geometry.get_z(xprime, yprime, zprime, pars.incl)

        rcyls[i] = rcyl
        zs[i] = z

        Deltav = v0 * 1e5 - DiskTracer.geometry.get_vlos(xprime, yprime, zprime, pars)
        Deltavs[i] = Deltav

        if zprime > 0
            tildeDeltav = DiskTracer.geometry.get_Delta_v_approx(xprime, yprime, zprime, pars, v0, z0s[1])
        else
            tildeDeltav = DiskTracer.geometry.get_Delta_v_approx(xprime, yprime, zprime, pars, v0, z0s[2])
        end
        tildeDeltavs[i] = tildeDeltav * 1e-5 #km/s

        DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)
        DeltaVs[i] = DeltaV * 1e-5 # km/s

        Deltas[i] = exp(-Deltav^2/DeltaV^2)

        # S = DiskTracer.model.Sigma(rcyl, pars)
        # H = (rcyl)^(0.5 * (3 - pars.q)) * ((10 * AU)^pars.q * kB * pars.T_10 / (mu_gas * m_H * G * pars.M_star * M_sun))^(1/2)
        # rho = S/(sqrt(2. * pi) * H) * exp(-0.5 * (z/H)^2)
        #
        # rhos[i] = rho # DiskTracer.model.rho(rcyl, z, pars)

        # Ss[i] = DiskTracer.model.S_nu(rcyl, z, mol.nu_0, pars)

        Us[i] = DiskTracer.model.Upsilon_nu(rcyl, z, mol.nu_0, pars, mol)
        tildeUs[i] = Upsilon_amp * exp(-0.5 * (zprime - zprime_midplane)^2/Upsilon_sigma^2)

        # println(Us[i], " ", tildeUs[i])

        alphas[i] = Us[i] * exp(-Deltav^2/DeltaV^2)
        tildealphas[i] = tildeUs[i] * exp(-tildeDeltavs[i]^2 / DeltaVs[i]^2)

        if zprime > 0
            tildealphasG[i] = a_alpha1 / sqrt(2pi * sigma_c2) * exp(-0.5 * (zprime - mu_c1)^2/sigma_c2)
        else
            tildealphasG[i] = a_alpha2 / sqrt(2pi * sigma_c2) * exp(-0.5 * (zprime - mu_c2)^2/sigma_c2)
        end
    end

    # get bounding zps
    zbounds = DiskTracer.geometry.get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(zbounds) == 2
        zstart, zend = zbounds
        println("H bound verified merged: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, zstart, zend, pars, 0.5))
    else
        z1start, z1end, z2start, z2end = zbounds
        println("H bound verified front-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z1start, z1end, pars, 0.5))
        println("H bound verified rear-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z2start, z2end, pars, 0.5))
    end

    # hbounds, derivs = DiskTracer.geometry.get_bounding_scale_heights(xprime, yprime, pars, 1.0)

    p1 = plot(zps./AU, rcyls./AU, ylabel=L"$r_\mathrm{cyl}$ [AU]")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p2 = plot(zps./AU, zs./AU, ylabel=L"$z$ [AU]")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p3 = plot(zps./AU, Deltavs .* 1e-5, ylabel=L"$\Delta v$")
    plot!(zps./AU, tildeDeltavs, ylims=(-1., 1))

    p4 = plot(zps./AU, Deltas, xlabel=L"$z^\prime$", ylabel=L"$\exp(-\Delta v^2/\Delta V^2)$")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))
    vline!((z0s./AU)', line=([:black :black]))

    p5 = plot(zps./AU, Us, xlabel=L"$z^\prime$", ylabel=L"$\Upsilon$") #, yaxis=(scale=:log10))

    p55 = plot(zps./AU, tildeUs, ylabel=L"$\tilde{\Upsilon}$") #, yaxis=(scale=:log10))
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    # p6 = plot(zps./AU, rhos, ylabel=L"$\rho$")

    # p45 = plot(zps./AU, DeltaVs, xlabel=L"$z^\prime$", ylabel=L"$\Delta V$")

    p7 = plot(zps./AU, alphas, xlabel=L"$z^\prime$", ylabel=L"$\alpha$")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))
    plot!(zps./AU, tildealphas, ylabel=L"$\tilde{\alpha}$")
    plot!(zps./AU, tildealphasG, ylabel=L"$\tilde{\alpha}$")

    p8 = plot(zps./AU, log10.(alphas) - log10.(tildealphas), ylim=(-2,2)) # ylabel=L"$\alpha/\tilde{\alpha}$")
    p9 = plot(zps./AU, log10.(alphas) - log10.(tildealphasG), ylim=(-2,2))

    # println(extrema(alphas./tildealphas))

    p = plot(p1,p2,p3,p4,p5,p55,p7,p8,p9, layout=(9,1), legend=false, size=(800,1200))

    savefig(p, "approx_ray_plot.png")

end

plot_ray(100.*AU, -100.0*AU)


text:
module Mattermost

using JSON2

export OutgoingWebhookRequest, OutgoingWebhookResponse, OutgoingWebhookResponseType

"""
    OutgoingWebhookRequest
Request type for outgoing webhooks. This is just a subset of the attributes.
See https://developers.mattermost.com/integrate/outgoing-webhooks/ for a
somewhat complete description of the JSON format.
"""
struct OutgoingWebhookRequest
    text::String
    user_name::String
    OutgoingWebhookRequest(; text, user_name, kwargs...) = new(text, user_name)
end

# format specification for OutgoingWebhookRequest
# keywordargs: parse arbitrarily ordered and extra fields
JSON2.@format OutgoingWebhookRequest keywordargs begin
end

"""
    OutgoingWebhookResponseType

Enum for `response_type`. Can be post (default if not specified) or comment.
"""
@enum OutgoingWebhookResponseType post comment

"""
    OutgoingWebhookResponse
Response type for outgoing webhooks. This is just a subset of the attributes.
See https://developers.mattermost.com/integrate/outgoing-webhooks/ for a
somewhat complete description of the JSON format.
"""
struct OutgoingWebhookResponse
    text::String
    response_type::Union{OutgoingWebhookResponseType, Nothing}

    OutgoingWebhookResponse(text, response_type = nothing) = new(text, response_type)
end

# format specification for OutgoingWebhookResponse
JSON2.@format OutgoingWebhookResponse begin
    # omitempty=true: don't write JSON null value for nothing
    response_type => (omitempty = true,)
end

end


text:
using PyCall

run(`$(PyCall.pyprogramname) -m pip install git+https://github.com/rgerum/pylustrator@61927a7be1ca2ea03f9378f85432c69e302f2c34`)


text:
## TODO: clean up this file

include("../solvers/tridiagonalsolver.jl")
solver(::TridiagonalProlateStyle, ap, A; scaling = Zt_scaling_factor(dictionary(ap), A)
    , options...) =
    FE_TridiagonalSolver(A, scaling; options...)
# They have to do with a normalization of the
# sampling operators.
Zt_scaling_factor(S::Dictionary, A) = length(supergrid(grid(dest(A))))
Zt_scaling_factor(S::DerivedDict, A) = Zt_scaling_factor(superdict(S), A)
Zt_scaling_factor(S::ChebyshevT, A) = length(supergrid(grid(dest(A))))/2

oversampling_grid(dict::ExtensionFrame, L) = subgrid(oversampling_grid(superdict(dict),L), support(dict))
# BasisFunctions.uniformweights(grid::SubGrid) =

discretemeasure(ss::SamplingStyle, platform::ExtensionFramePlatform, param, ap_old; options...) =
    restrict(discretemeasure(ss, platform.basisplatform, param,
        (ap=approximationproblem(platform.basisplatform, param); cache!(ap,samplingparameter, samplingparameter(ap_old));ap)
            ; options...), platform.domain)

azdual(sstyle::SamplingStyle, platform::ExtensionFramePlatform, param, L, measure::Measure; options...) =
   extensionframe(azdual(sstyle, platform.basisplatform, param, L, supermeasure(measure); options...), platform.domain)

correct_sampling_parameter(platform::ExtensionFramePlatform, param, L_trial; options...) =
   correct_sampling_parameter(platform.basisplatform, param, L_trial; options...)

azdual(sstyle::SamplingStyle, platform::AugmentationPlatform, param, L, measure::Measure; options...) =
   MultiDict([azdual(sstyle, platform.basis, param, L, measure; options...), platform.functions])

oversampling_grid(samplingstyle::SamplingStyle, platform::AugmentationPlatform, param, L; options...) =
   oversampling_grid(samplingstyle, platform.basis, param, L; options...)

correct_sampling_parameter(platform::AugmentationPlatform, param, L; options...) =
    correct_sampling_parameter(platform.basis, param, L; options...)

function azdual(sstyle::SamplingStyle, platform::WeightedSumPlatform, param, L, measure::Measure; options...)
    denom = (x...)->sum(map(w->abs(w(x...))^2, platform.weights))
    # TODO: discuss, what is the relation between param, L of a platform and platform.P
    MultiDict([((x...)->(platform.weights[j](x...)/denom(x...))) * azdual(sstyle, platform.P, param, L, measure; options...) for j=1:length(platform.weights)])
end

oversampling_grid(samplingstyle::SamplingStyle, platform::WeightedSumPlatform, param::NTuple, L; options...) =
  oversampling_grid(samplingstyle, platform.P, param[1], L; options...)

azdual(sstyle::SamplingStyle, dict::ExtensionFrame, L, measure::Measure; options...) =
   extensionframe(support(dict), azdual(sstyle, superdict(dict), L, supermeasure(measure); options...),)

function azdual(dict::MultiDict, measure::Measure; options...)
    dictionary = dict.dicts[1].superdict
    weights = map(weightfunction, components(dict))
    denom = (x...)->sum(map(w->abs(w(x...))^2, weights))
    MultiDict([((x...)->(weights[j](x...)/denom(x...))) * azdual(dictionary, measure; options...) for j=1:length(weights)])
end

deduce_samplingparameter(ss::OversamplingStyle, platform::WeightedSumPlatform, param::NTuple; options...) =
    mix_samplingparameters(platform, map(parami->deduce_oversampling_parameter(ss, platform.P[parami]; options...), param))

mix_samplingparameters(platform::WeightedSumPlatform, Ls::NTuple) =
    round.(Int,(sum(prod.(Ls))/prod(Ls[1]))^(1/length(Ls[1]))     .*Ls[1])

correct_sampling_parameter(platform::WeightedSumPlatform, param, L; options...) =
    correct_sampling_parameter(platform.P, param, L; options...)


text:
abstract type Integrator end

Base.broadcastable(x::Integrator) = Ref(x)

struct NoIntegrator <: Integrator end

struct VVIntegrator <: Integrator
    hdtsqr::Float64
    function VVIntegrator(dt::Float64)
        hdtsqr=dt*dt/2.0
        new(hdtsqr::Float64)
    end
end # Velocity Verlet


text:
#=----------------------------------------------------------------------------

Copyright (c) 2015 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

----------------------------------------------------------------------------=#

import Images
export relief

#---------------------------------------------------------------------
"""
relief -  Generates relief shaded image

```
Usage:  shadeimg = relief(img, azimuth, elevation, gradscale, rgbimg)

Arguments: img - Image/heightmap to be relief shaded. ::Array{T<:Real,2}
       azimuth - Of light direction in degrees. Zero azimuth points
                 upwards and increases clockwise. Defaults to 45.
     elevation - Of light direction in degrees. Defaults to 45.
     gradscale - Scaling to apply to the surface gradients.  If the shading
                 is excessive decrease the scaling. Try successive doubling
                 or halving to find a good value.
        rgbimg - Optional RGB image to which the shading pattern derived
                 from 'img' is applied.   Alternatively, rgbimg can be a
                 colour map of type ::Array{ColorTypes.RGBA{Float64},1}
                 obtained from cmap().  This colour map is applied to the input
                 image/heightmap in order to obtain a RGB image to which
                 the shading pattern is applied.
```

Lambertian shading is used to form the relief image.  This obtained from the
cosine of the angle between the surface normal and light direction.  Note that
shadows are ignored.  Thus a small feature that might otherwise be in the
shadow of a nearby large structure is rendered as if the large feature was not
there.

See also: cmap, applycolourmap
"""
function relief(img::Array{T1,2}, az::Real=45, el::Real=45,
                gradscale::Real=1, rgbimg::Array{T2,3}=Array{Float64}(undef, 0,0,0)) where {T1<:Real,T2<:Real}
    # April    2014 Original MATLAB version
    # December 2015 Ported to Julia

    (rows, cols) = size(img)

    if !isempty(rgbimg)
        (r,c,chan) = size(rgbimg)
        if rows != r || cols != c || chan != 3
            error("rgbimg must be the same size as img")
        end
    end

    # Obtain surface normals of img
    loggrad = "lin"
    (n1, n2, n3) = surfacenormals(img, gradscale, loggrad)

    # Convert azimuth and elevation to a lighting direction vector.  Note that
    # the vector is constructed so that an azimuth of 0 points upwards and
    # increases clockwise.
    az = az/180*pi
    el = el/180*pi
    I = [cos(el)*sin(az), cos(el)*cos(az), sin(el)]
    I = I./norm(I)  # Ensure I is a unit vector

    # Generate Lambertian shading via the dot product between surface normals
    # and the light direction.  Note that the product with n2 is negated to
    # account for the image +ve y increasing downwards.
    shading = I[1]*n1 - I[2]*n2 + I[3]*n3

    # Remove -ve shading values which are generated by surface normals pointing
    # away from the light source.
    shading[shading .< 0] .= 0

    # If no RGB image has been supplied just return the raw shading image
    if isempty(rgbimg)
        shadeimg = shading

    else  #  Apply shading to the RGB image supplied
        shadeimg = zeros(size(rgbimg))

        for n = 1:3
            shadeimg[:,:,n]  = rgbimg[:,:,n].*shading
        end
    end

    # ** Resolve issue with RGB image being double or uint8

    return shadeimg
end

# Case where a colour map is supplied rather than a RGB img
function relief(img::Array{T1,2}, az::Real, el::Real,
                gradscale::Real,
                rgbmap::Array{ColorTypes.RGB{Float64},1}) where T1<:Real

    # Apply colour map to input image and use this as rgbimg
    rgbimg = applycolourmap(img, rgbmap)
    return relief(img, az, el, gradscale, rgbimg)
end

# Case where a colour map is supplied rather than a RGB img
function relief(img::Array{T1,2}, az::Real, el::Real,
                gradscale::Real,
                rgbmap::Array{ColorTypes.RGBA{Float64},1}) where T1<:Real

    # Apply colour map to input image and use this as rgbimg
    rgbimg = applycolourmap(img, rgbmap)
    return relief(img, az, el, gradscale, rgbimg)
end

#---------------------------------------------------------------------------
# Compute image/heightmap surface normals
#
# (n1, n2, n3) = surfacenormals(img, gradscale, loggrad)

function surfacenormals(img::Array{T,2}, gradscale::Real, loggrad::String) where T<:Real

    loggrad = lowercase(loggrad)

    # Compute partial derivatives of z.
    # p = dz/dx, q = dz/dy
    (q,p) = Images.imgradients(img, KernelFactors.ando3)
    p = p*gradscale
    q = q*gradscale

    # If specified take logs of gradient.
    # Note that taking the log of the surface gradients will produce a surface
    # that is not integrable (probably only of theoretical interest)
    if loggrad == "log"
        p = sign(p).*log1p(abs(p))
        q = sign(q).*log1p(abs(q))
    elseif loggrad == "loglog"
        p = sign(p).*log1p(log1p(abs(p)))
        q = sign(q).*log1p(log1p(abs(q)))
    elseif loggrad == "logloglog"
        p = sign(p).*log1p(log1p(log1p(abs(p))))
        q = sign(q).*log1p(log1p(log1p(abs(q))))
    end

    # Generate surface unit normal vectors. Components stored in n1, n2
    # and n3
    mag = sqrt.(1 .+ p.^2 .+ q.^2)
    n1 = -p./mag
    n2 = -q./mag
    n3 =  1 ./ mag

    return n1, n2, n3
end


text:
import PhysicsEngine

using HTTP

@testset "Line Breaking" begin
  request = HTTP.request("GET", "http://www.unicode.org/Public/UCD/latest/ucd/auxiliary/LineBreakTest.txt")
  body = String(request.body)
  working = true
  line_n = 1
  fails = 0
  total = 0

  # Bad tests to skip
  skip_lines = [
    1141,
    1145,
    1309,
    1313,
    2981,
    4497,
    4665,
    5165
  ]

  for line in split(body, '\n')
    main_string = split(line, '#')[1]
    tokens = split(main_string)
    breaks = UInt32[]
    seq = ""
    skip = false
    for token in tokens
      if token == "×"
        push!(breaks, UInt32(1))
      elseif token == "÷"
        push!(breaks, UInt32(0))
      else
        if parse(UInt32, token, base=16) > 128
          skip = true
        end
        seq = seq * string(Char(parse(UInt32, token, base=16)))
      end
    end

    if line_n in skip_lines
      skip = true
    end

    if skip
      line_n += 1
      continue
    end

    if length(tokens) > 0
      total += 1
    end

    generatedBreaksPtr = ccall((:lineBreaks, PhysicsEngine.@fullLibraryPath), Ptr{UInt32}, (Cstring,), seq)
    generatedBreaks = UInt32[]
    # Make an array
    for i = 1:(length(seq) + 1)
      action = Base.unsafe_load(generatedBreaksPtr, i)
      if action == 2
        action = 0
      end
      push!(generatedBreaks, UInt32(action))
    end

    broken = false
    for (expected, actual) in zip(breaks, generatedBreaks)
      if expected != actual
        broken = true
        working = false
      end
    end

    if broken
      fails += 1
      print("$line_n - ")
      index = 1
      for c in seq
        ex = breaks[index]
        ac = generatedBreaks[index]
        print(" [$ex|$ac] ")
        print(UInt32(c))
        #print("('$c')")
        print(" ")
        index += 1
      end
      ex = breaks[index]
      ac = generatedBreaks[index]
      print(" [$ex|$ac]")
      print("\n")
    end
    line_n += 1
  end
  print("Fails: $fails / $total \n")
  @test working == true
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
module Measures

export Measure, Length, AbsoluteLength, BoundingBox, AbsoluteBox, Absolute2DBox,
       Absolute3DBox, Vec, Vec2, Vec3, AbsoluteVec, AbsoluteVec2, AbsoluteVec3,
       isabsolute, resolve,
        mm, cm, inch, pt, width, height

abstract type Measure end

include("operations.jl")
include("length.jl")
include("point.jl")
include("boundingbox.jl")
include("resolution.jl")

Base.show(out::IO, x::Length{U}) where U = print(out, x.value, U)
Base.show(out::IO, x::Neg) = print(out, "-", x.a)
Base.show(out::IO, x::Add) = print(out, x.a, " + ", x.b)
Base.show(out::IO, x::Min) = print(out, "min(", x.a, ", ", x.b, ")")
Base.show(out::IO, x::Max) = print(out, "max(", x.a, ", ", x.b, ")")
Base.show(out::IO, x::Div) = print(out,  x.a, " / ", x.b)
Base.show(out::IO, x::Mul) = print(out,  x.a, " * ", x.b)

Broadcast.broadcastable(x::T) where T<:Measure = Ref(x)

end # module Measures


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function DataFile(TSTART,TSTOP,Ts)
# ----------------------------------------------------------------------------------- #
# DataFile.jl was generated using the Kwatee code generation system.
# DataFile: Stores model parameters as key - value pairs in a Julia Dict()
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 04-19-2016 10:25:05
#
# Input arguments:
# TSTART  - Time start
# TSTOP  - Time stop
# Ts - Time step
#
# Return arguments:
# data_dictionary  - Data dictionary instance (holds model parameters)
# ----------------------------------------------------------------------------------- #

# Load the stoichiometric matrix -
S = float(open(readdlm,"/Users/dbassen/work/KwateeServer-v1.0/network/Network.dat"));
(NSPECIES,NREACTIONS) = size(S);

# How many genes do we have in the model? -
number_of_genes = 22;

# Formulate the initial condition array -
initial_condition_array = Float64[];
push!(initial_condition_array,1.0);	#	1	gene_Smad
push!(initial_condition_array,1.0);	#	2	gene_Ecadherin
push!(initial_condition_array,1.0);	#	3	gene_Bcatenin
push!(initial_condition_array,1.0);	#	4	gene_PLCy
push!(initial_condition_array,1.0);	#	5	gene_PKC
push!(initial_condition_array,1.0);	#	6	gene_RAF
push!(initial_condition_array,1.0);	#	7	gene_RAS
push!(initial_condition_array,1.0);	#	8	gene_AP1_SP1
push!(initial_condition_array,1.0);	#	9	gene_tubulin
push!(initial_condition_array,1.0);	#	10	gene_SNAIL_SLUG
push!(initial_condition_array,1.0);	#	11	gene_Calcineurin
push!(initial_condition_array,1.0);	#	12	gene_GSK3_APC_AXIN
push!(initial_condition_array,1.0);	#	13	gene_PI3K
push!(initial_condition_array,0.0);	#	14	mRNA_Smad
push!(initial_condition_array,0.0);	#	15	mRNA_Ecadherin
push!(initial_condition_array,0.0);	#	16	mRNA_Bcatenin
push!(initial_condition_array,0.0);	#	17	mRNA_PLCy
push!(initial_condition_array,0.0);	#	18	mRNA_PKC
push!(initial_condition_array,0.0);	#	19	mRNA_RAF
push!(initial_condition_array,0.0);	#	20	mRNA_RAS
push!(initial_condition_array,0.0);	#	21	mRNA_AP1_SP1
push!(initial_condition_array,0.0);	#	22	mRNA_Calcineurin
push!(initial_condition_array,0.0);	#	23	mRNA_SNAIL_SLUG
push!(initial_condition_array,0.0);	#	24	mRNA_tubulin
push!(initial_condition_array,0.0);	#	25	mRNA_GSK3_APC_AXIN
push!(initial_condition_array,0.0);	#	26	mRNA_PI3K
push!(initial_condition_array,0.0);	#	27	PLCy
push!(initial_condition_array,0.0);	#	28	Active_PLCy
push!(initial_condition_array,0.0);	#	29	Calcineurin
push!(initial_condition_array,0.0);	#	30	Active_Calcineurin
push!(initial_condition_array,0.0);	#	31	NFATc
push!(initial_condition_array,0.0);	#	32	Active_NFATc
push!(initial_condition_array,0.0);	#	33	PKC
push!(initial_condition_array,0.0);	#	34	Active_PKC
push!(initial_condition_array,0.0);	#	35	RAF
push!(initial_condition_array,0.0);	#	36	RAF_P
push!(initial_condition_array,0.0);	#	37	Smad
push!(initial_condition_array,0.0);	#	38	pSmad
push!(initial_condition_array,0.0);	#	39	PI3K
push!(initial_condition_array,0.0);	#	40	PI3K_P
push!(initial_condition_array,0.0);	#	41	GSK3_APC_AXIN
push!(initial_condition_array,0.0);	#	42	GSK3_APC_AXIN_P
push!(initial_condition_array,0.0);	#	43	RAS
push!(initial_condition_array,0.0);	#	44	RAS_GTP
push!(initial_condition_array,0.0);	#	45	MAPK
push!(initial_condition_array,0.0);	#	46	pMAPK
push!(initial_condition_array,0.0);	#	47	AP1_SP1
push!(initial_condition_array,0.0);	#	48	AP1_SP1_P
push!(initial_condition_array,0.0);	#	49	Ecadherin
push!(initial_condition_array,0.0);	#	50	Bcatenin
push!(initial_condition_array,0.0);	#	51	Ecadherin_Bcatenin
push!(initial_condition_array,0.0);	#	52	GSK3_APC_AXIN_Bcatenin
push!(initial_condition_array,0.0);	#	53	TCF4
push!(initial_condition_array,0.0);	#	54	Active_TCF4
push!(initial_condition_array,0.0);	#	55	LEF1
push!(initial_condition_array,0.0);	#	56	Active_LEF1
push!(initial_condition_array,0.0);	#	57	PPase
push!(initial_condition_array,0.0);	#	58	YREG1
push!(initial_condition_array,0.0);	#	59	SNAIL_SLUG
push!(initial_condition_array,0.0);	#	60	vimentin
push!(initial_condition_array,0.0);	#	61	tubulin
push!(initial_condition_array,0.0);	#	62	RNAP
push!(initial_condition_array,0.0);	#	63	RIBOSOME
push!(initial_condition_array,0.0);	#	64	U0126
push!(initial_condition_array,0.0);	#	65	TGFB3_Ab
push!(initial_condition_array,0.0);	#	66	ss_virus
push!(initial_condition_array,0.0);	#	67	smad_virus
push!(initial_condition_array,0.0);	#	68	lef1_virus
push!(initial_condition_array,1.0);	#	69	gene_TCF4
push!(initial_condition_array,1.0);	#	70	gene_TGFB3
push!(initial_condition_array,1.0);	#	71	gene_YREG1
push!(initial_condition_array,1.0);	#	72	gene_LEF1
push!(initial_condition_array,1.0);	#	73	gene_vimentin
push!(initial_condition_array,1.0);	#	74	gene_R
push!(initial_condition_array,1.0);	#	75	gene_VEGFR
push!(initial_condition_array,1.0);	#	76	gene_NFATc
push!(initial_condition_array,1.0);	#	77	gene_MAPK
push!(initial_condition_array,0.0);	#	78	mRNA_TCF4
push!(initial_condition_array,0.0);	#	79	mRNA_TGFB3
push!(initial_condition_array,0.0);	#	80	mRNA_YREG1
push!(initial_condition_array,0.0);	#	81	mRNA_LEF1
push!(initial_condition_array,0.0);	#	82	mRNA_vimentin
push!(initial_condition_array,0.0);	#	83	mRNA_R
push!(initial_condition_array,0.0);	#	84	mRNA_VEGFR
push!(initial_condition_array,0.0);	#	85	mRNA_NFATc
push!(initial_condition_array,0.0);	#	86	mRNA_MAPK
push!(initial_condition_array,0.0);	#	87	R
push!(initial_condition_array,0.0);	#	88	TGFB12
push!(initial_condition_array,0.0);	#	89	TGFB3
push!(initial_condition_array,0.0);	#	90	VEGF
push!(initial_condition_array,0.0);	#	91	TGFB12_R
push!(initial_condition_array,0.0);	#	92	TGFB3_R
push!(initial_condition_array,0.0);	#	93	VEGFR
push!(initial_condition_array,0.0);	#	94	VEGF_VEGFR

# Formulate the time constant array -
time_constant_array = Float64[];
push!(time_constant_array,1.0);	#	1	 time constant: gene_Smad
push!(time_constant_array,1.0);	#	2	 time constant: gene_Ecadherin
push!(time_constant_array,1.0);	#	3	 time constant: gene_Bcatenin
push!(time_constant_array,1.0);	#	4	 time constant: gene_PLCy
push!(time_constant_array,1.0);	#	5	 time constant: gene_PKC
push!(time_constant_array,1.0);	#	6	 time constant: gene_RAF
push!(time_constant_array,1.0);	#	7	 time constant: gene_RAS
push!(time_constant_array,1.0);	#	8	 time constant: gene_AP1_SP1
push!(time_constant_array,1.0);	#	9	 time constant: gene_tubulin
push!(time_constant_array,1.0);	#	10	 time constant: gene_SNAIL_SLUG
push!(time_constant_array,1.0);	#	11	 time constant: gene_Calcineurin
push!(time_constant_array,1.0);	#	12	 time constant: gene_GSK3_APC_AXIN
push!(time_constant_array,1.0);	#	13	 time constant: gene_PI3K
push!(time_constant_array,0.1);	#	14	 time constant: mRNA_Smad
push!(time_constant_array,0.1);	#	15	 time constant: mRNA_Ecadherin
push!(time_constant_array,0.1);	#	16	 time constant: mRNA_Bcatenin
push!(time_constant_array,0.1);	#	17	 time constant: mRNA_PLCy
push!(time_constant_array,0.1);	#	18	 time constant: mRNA_PKC
push!(time_constant_array,0.1);	#	19	 time constant: mRNA_RAF
push!(time_constant_array,0.1);	#	20	 time constant: mRNA_RAS
push!(time_constant_array,0.1);	#	21	 time constant: mRNA_AP1_SP1
push!(time_constant_array,0.1);	#	22	 time constant: mRNA_Calcineurin
push!(time_constant_array,0.1);	#	23	 time constant: mRNA_SNAIL_SLUG
push!(time_constant_array,0.1);	#	24	 time constant: mRNA_tubulin
push!(time_constant_array,0.1);	#	25	 time constant: mRNA_GSK3_APC_AXIN
push!(time_constant_array,0.1);	#	26	 time constant: mRNA_PI3K
push!(time_constant_array,1.0);	#	27	 time constant: PLCy
push!(time_constant_array,1.0);	#	28	 time constant: Active_PLCy
push!(time_constant_array,1.0);	#	29	 time constant: Calcineurin
push!(time_constant_array,1.0);	#	30	 time constant: Active_Calcineurin
push!(time_constant_array,1.0);	#	31	 time constant: NFATc
push!(time_constant_array,1.0);	#	32	 time constant: Active_NFATc
push!(time_constant_array,1.0);	#	33	 time constant: PKC
push!(time_constant_array,1.0);	#	34	 time constant: Active_PKC
push!(time_constant_array,1.0);	#	35	 time constant: RAF
push!(time_constant_array,1.0);	#	36	 time constant: RAF_P
push!(time_constant_array,1.0);	#	37	 time constant: Smad
push!(time_constant_array,1.0);	#	38	 time constant: pSmad
push!(time_constant_array,1.0);	#	39	 time constant: PI3K
push!(time_constant_array,1.0);	#	40	 time constant: PI3K_P
push!(time_constant_array,1.0);	#	41	 time constant: GSK3_APC_AXIN
push!(time_constant_array,1.0);	#	42	 time constant: GSK3_APC_AXIN_P
push!(time_constant_array,1.0);	#	43	 time constant: RAS
push!(time_constant_array,1.0);	#	44	 time constant: RAS_GTP
push!(time_constant_array,1.0);	#	45	 time constant: MAPK
push!(time_constant_array,1.0);	#	46	 time constant: pMAPK
push!(time_constant_array,1.0);	#	47	 time constant: AP1_SP1
push!(time_constant_array,1.0);	#	48	 time constant: AP1_SP1_P
push!(time_constant_array,1.0);	#	49	 time constant: Ecadherin
push!(time_constant_array,1.0);	#	50	 time constant: Bcatenin
push!(time_constant_array,1.0);	#	51	 time constant: Ecadherin_Bcatenin
push!(time_constant_array,1.0);	#	52	 time constant: GSK3_APC_AXIN_Bcatenin
push!(time_constant_array,1.0);	#	53	 time constant: TCF4
push!(time_constant_array,1.0);	#	54	 time constant: Active_TCF4
push!(time_constant_array,1.0);	#	55	 time constant: LEF1
push!(time_constant_array,1.0);	#	56	 time constant: Active_LEF1
push!(time_constant_array,1.0);	#	57	 time constant: PPase
push!(time_constant_array,1.0);	#	58	 time constant: YREG1
push!(time_constant_array,1.0);	#	59	 time constant: SNAIL_SLUG
push!(time_constant_array,1.0);	#	60	 time constant: vimentin
push!(time_constant_array,1.0);	#	61	 time constant: tubulin
push!(time_constant_array,1.0);	#	62	 time constant: RNAP
push!(time_constant_array,1.0);	#	63	 time constant: RIBOSOME
push!(time_constant_array,1.0);	#	64	 time constant: U0126
push!(time_constant_array,1.0);	#	65	 time constant: TGFB3_Ab
push!(time_constant_array,1.0);	#	66	 time constant: ss_virus
push!(time_constant_array,1.0);	#	67	 time constant: smad_virus
push!(time_constant_array,1.0);	#	68	 time constant: lef1_virus
push!(time_constant_array,1.0);	#	69	 time constant: gene_TCF4
push!(time_constant_array,1.0);	#	70	 time constant: gene_TGFB3
push!(time_constant_array,1.0);	#	71	 time constant: gene_YREG1
push!(time_constant_array,1.0);	#	72	 time constant: gene_LEF1
push!(time_constant_array,1.0);	#	73	 time constant: gene_vimentin
push!(time_constant_array,1.0);	#	74	 time constant: gene_R
push!(time_constant_array,1.0);	#	75	 time constant: gene_VEGFR
push!(time_constant_array,1.0);	#	76	 time constant: gene_NFATc
push!(time_constant_array,1.0);	#	77	 time constant: gene_MAPK
push!(time_constant_array,0.1);	#	78	 time constant: mRNA_TCF4
push!(time_constant_array,0.1);	#	79	 time constant: mRNA_TGFB3
push!(time_constant_array,0.1);	#	80	 time constant: mRNA_YREG1
push!(time_constant_array,0.1);	#	81	 time constant: mRNA_LEF1
push!(time_constant_array,0.1);	#	82	 time constant: mRNA_vimentin
push!(time_constant_array,0.1);	#	83	 time constant: mRNA_R
push!(time_constant_array,0.1);	#	84	 time constant: mRNA_VEGFR
push!(time_constant_array,0.1);	#	85	 time constant: mRNA_NFATc
push!(time_constant_array,0.1);	#	86	 time constant: mRNA_MAPK
push!(time_constant_array,1.0);	#	87	 time constant: R
push!(time_constant_array,1.0);	#	88	 time constant: TGFB12
push!(time_constant_array,1.0);	#	89	 time constant: TGFB3
push!(time_constant_array,1.0);	#	90	 time constant: VEGF
push!(time_constant_array,1.0);	#	91	 time constant: TGFB12_R
push!(time_constant_array,1.0);	#	92	 time constant: TGFB3_R
push!(time_constant_array,1.0);	#	93	 time constant: VEGFR
push!(time_constant_array,1.0);	#	94	 time constant: VEGF_VEGFR

# Formulate the rate constant array -
rate_constant_array = Float64[];
push!(rate_constant_array,1.0);	# 1	TGFB12_binding: TGFB12+R =([])=> TGFB12_R
push!(rate_constant_array,1.0);	# 2	TGFB12_binding: TGFB12_R =([])=> TGFB12+R (reverse)
push!(rate_constant_array,1.0);	# 3	TGFB3_binding: TGFB3+R =([])=> TGFB3_R
push!(rate_constant_array,1.0);	# 4	TGFB3_binding: TGFB3_R =([])=> TGFB3+R (reverse)
push!(rate_constant_array,1.0);	# 5	VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR
push!(rate_constant_array,1.0);	# 6	VEGF_binding: VEGF_VEGFR =([])=> VEGF+VEGFR (reverse)
push!(rate_constant_array,1.0);	# 7	VEGF_VEGFR_PLCy: PLCy =(VEGF_VEGFR)=> Active_PLCy
push!(rate_constant_array,1.0);	# 8	PLCy_activation_Calcineurin: Calcineurin =(Active_PLCy)=> Active_Calcineurin
push!(rate_constant_array,1.0);	# 9	Calcineurin_activation_NFATc: NFATc =(Active_Calcineurin)=> Active_NFATc
push!(rate_constant_array,1.0);	# 10	PLCy_activatation_PKC: PKC =(Active_PLCy)=> Active_PKC
push!(rate_constant_array,1.0);	# 11	PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P
push!(rate_constant_array,1.0);	# 12	TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad
push!(rate_constant_array,1.0);	# 13	TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad
push!(rate_constant_array,1.0);	# 14	TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P
push!(rate_constant_array,1.0);	# 15	PI3K_phosphorylation_GSK3: GSK3_APC_AXIN =(PI3K_P)=> GSK3_APC_AXIN_P
push!(rate_constant_array,1.0);	# 16	TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
push!(rate_constant_array,1.0);	# 17	RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
push!(rate_constant_array,1.0);	# 18	RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
push!(rate_constant_array,1.0);	# 19	MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad
push!(rate_constant_array,1.0);	# 20	MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
push!(rate_constant_array,1.0);	# 21	Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
push!(rate_constant_array,1.0);	# 22	Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
push!(rate_constant_array,1.0);	# 23	Binding_GSK3_APC_AXIN: GSK3_APC_AXIN+Bcatenin =([])=> GSK3_APC_AXIN_Bcatenin
push!(rate_constant_array,1.0);	# 24	Binding_GSK3_APC_AXIN: GSK3_APC_AXIN_Bcatenin =([])=> GSK3_APC_AXIN+Bcatenin (reverse)
push!(rate_constant_array,1.0);	# 25	TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
push!(rate_constant_array,1.0);	# 26	TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
push!(rate_constant_array,1.0);	# 27	LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
push!(rate_constant_array,1.0);	# 28	LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)
push!(rate_constant_array,1.0);	# 29	deactivation_pSmad: pSmad =(PPase)=> Smad
push!(rate_constant_array,1.0);	# 30	deactivation_RAF_P: RAF_P =(PPase)=> RAF
push!(rate_constant_array,1.0);	# 31	deactivation_RAS_GTP: RAS_GTP =([])=> RAS
push!(rate_constant_array,1.0);	# 32	deactivation_pMAPK: pMAPK =(PPase)=> MAPK
push!(rate_constant_array,1.0);	# 33	deactivation_AP1_SP1: AP1_SP1_P =(PPase)=> AP1_SP1
push!(rate_constant_array,1.0);	# 34	deactivation_Active_PLCy: Active_PLCy =([])=> PLCy
push!(rate_constant_array,1.0);	# 35	deactivation_Active_PKC: Active_PKC =([])=> PKC
push!(rate_constant_array,1.0);	# 36	deactivation_Active_Calcineurin: Active_Calcineurin =([])=> Calcineurin
push!(rate_constant_array,1.0);	# 37	deactivation_Active_NFATc: Active_NFATc =([])=> NFATc
push!(rate_constant_array,1.0);	# 38	deactivation_GSK3_APC_AXIN_P: GSK3_APC_AXIN_P =(PPase)=> GSK3_APC_AXIN
push!(rate_constant_array,0.1);	# 39	degradation_TGFB12_R: TGFB12_R =([])=> []
push!(rate_constant_array,0.1);	# 40	degradation_TGFB3_R: TGFB3_R =([])=> []
push!(rate_constant_array,0.1);	# 41	degradation_TGFB3: TGFB3 =([])=> []
push!(rate_constant_array,0.1);	# 42	Bcatenin_degradation: Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 43	Ecadherin_degrdation: Ecadherin =([])=> []
push!(rate_constant_array,0.1);	# 44	Smad_degradation: Smad =([])=> []
push!(rate_constant_array,0.1);	# 45	TCF4_degradation: TCF4 =([])=> []
push!(rate_constant_array,0.1);	# 46	Active_TCF4_degradation: Active_TCF4 =([])=> []
push!(rate_constant_array,0.1);	# 47	YREG1_degradation: YREG1 =([])=> []
push!(rate_constant_array,0.1);	# 48	SNAIL_SLUG_degradation: SNAIL_SLUG =([])=> []
push!(rate_constant_array,0.1);	# 49	LEF1_degradation: LEF1 =([])=> []
push!(rate_constant_array,0.1);	# 50	Active_LEF1_degradation: Active_LEF1 =([])=> []
push!(rate_constant_array,0.1);	# 51	vimentin_degradation: vimentin =([])=> []
push!(rate_constant_array,0.1);	# 52	VEGFR_degradation: VEGFR =([])=> []
push!(rate_constant_array,0.1);	# 53	NFATc_degradation: NFATc =([])=> []
push!(rate_constant_array,0.1);	# 54	Active_NFATc_degradation: Active_NFATc =([])=> []
push!(rate_constant_array,0.1);	# 55	degradation_MAPK: MAPK =([])=> []
push!(rate_constant_array,0.1);	# 56	degradation_R: R =([])=> []
push!(rate_constant_array,0.1);	# 57	Ecadherin_Bcatenin_degradation: Ecadherin_Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 58	AP1_SP1_degradation: AP1_SP1 =([])=> []
push!(rate_constant_array,0.1);	# 59	AP1_SP1_P_degradation: AP1_SP1_P =([])=> []
push!(rate_constant_array,0.1);	# 60	RAS_degradation: RAS =([])=> []
push!(rate_constant_array,0.1);	# 61	RAF_degradation: RAF =([])=> []
push!(rate_constant_array,0.1);	# 62	RAS_GTP_degradation: RAS_GTP =([])=> []
push!(rate_constant_array,0.1);	# 63	RAF_P_degradation: RAF_P =([])=> []
push!(rate_constant_array,0.1);	# 64	pMAPK_degradation: pMAPK =([])=> []
push!(rate_constant_array,0.1);	# 65	PLCy_degradation: PLCy =([])=> []
push!(rate_constant_array,0.1);	# 66	Active_PLCy_degradation: Active_PLCy =([])=> []
push!(rate_constant_array,0.1);	# 67	PKC_degradation: PKC =([])=> []
push!(rate_constant_array,0.1);	# 68	Active_PKC_degradation: Active_PKC =([])=> []
push!(rate_constant_array,0.1);	# 69	Calcineurin_degradation: Calcineurin =([])=> []
push!(rate_constant_array,0.1);	# 70	Active_Calcineurin_degradation: Active_Calcineurin =([])=> []
push!(rate_constant_array,0.1);	# 71	degradation_GSK3_APC_AXIN_degradation: GSK3_APC_AXIN =([])=> []
push!(rate_constant_array,0.1);	# 72	degradation_GSK3_APC_AXIN_P_degradation: GSK3_APC_AXIN_P =([])=> []
push!(rate_constant_array,0.1);	# 73	degradation_GSK3_APC_AXIN_Bcatenin_degradation: GSK3_APC_AXIN_Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 74	tubulin_degradation: tubulin =([])=> []
push!(rate_constant_array,0.1);	# 75	mRNA_Smad_degradation: mRNA_Smad =([])=> []
push!(rate_constant_array,0.1);	# 76	mRNA_Ecadherin_degradation: mRNA_Ecadherin =([])=> []
push!(rate_constant_array,0.1);	# 77	mRNA_Bcatenin_degradation: mRNA_Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 78	mRNA_TCF4_degradation: mRNA_TCF4 =([])=> []
push!(rate_constant_array,0.1);	# 79	mRNA_TGFB3_degradation: mRNA_TGFB3 =([])=> []
push!(rate_constant_array,0.1);	# 80	mRNA_YREG1_degradation: mRNA_YREG1 =([])=> []
push!(rate_constant_array,0.1);	# 81	mRNA_LEF1_degradation: mRNA_LEF1 =([])=> []
push!(rate_constant_array,0.1);	# 82	mRNA_vimentin_degradation: mRNA_vimentin =([])=> []
push!(rate_constant_array,0.1);	# 83	mRNA_VEGFR_degradation: mRNA_VEGFR =([])=> []
push!(rate_constant_array,0.1);	# 84	mRNA_NFATc_degradation: mRNA_NFATc =([])=> []
push!(rate_constant_array,0.1);	# 85	mRNA_R_degradation: mRNA_R =([])=> []
push!(rate_constant_array,0.1);	# 86	mRNA_MAPK_degradation: mRNA_MAPK =([])=> []
push!(rate_constant_array,0.1);	# 87	mRNA_tubulin_degradation: mRNA_tubulin =([])=> []
push!(rate_constant_array,0.1);	# 88	mRNA_SNAIL_SLUG_degradation: mRNA_SNAIL_SLUG =([])=> []
push!(rate_constant_array,0.1);	# 89	mRNA_RAS_degradation: mRNA_RAS =([])=> []
push!(rate_constant_array,0.1);	# 90	mRNA_RAF_degradation: mRNA_RAF =([])=> []
push!(rate_constant_array,0.1);	# 91	mRNA_PLCy_degradation: mRNA_PLCy =([])=> []
push!(rate_constant_array,0.1);	# 92	mRNA_PKC_degradation: mRNA_PKC =([])=> []
push!(rate_constant_array,0.1);	# 93	mRNA_Calcineurin_degradation: mRNA_Calcineurin =([])=> []
push!(rate_constant_array,0.1);	# 94	mRNA_AP1_SP1_degradation: mRNA_AP1_SP1 =([])=> []
push!(rate_constant_array,0.1);	# 95	VEGF_VEGFR_degradation: VEGF_VEGFR =([])=> []
push!(rate_constant_array,0.1);	# 96	generation_RNAP: RNAP =([])=> []
push!(rate_constant_array,1.0);	# 97	generation_RNAP: [] =([])=> RNAP (reverse)
push!(rate_constant_array,0.1);	# 98	generation_RIBOSOME: RIBOSOME =([])=> []
push!(rate_constant_array,1.0);	# 99	generation_RIBOSOME: [] =([])=> RIBOSOME (reverse)
push!(rate_constant_array,0.1);	# 100	generation_PPase: PPase =([])=> []
push!(rate_constant_array,1.0);	# 101	generation_PPase: [] =([])=> PPase (reverse)
push!(rate_constant_array,0.1);	# 102	addition_U0126: U0126 =([])=> []
push!(rate_constant_array,1.0);	# 103	addition_U0126: [] =([])=> U0126 (reverse)
push!(rate_constant_array,0.1);	# 104	addition_TGFB3_Ab: TGFB3_Ab =([])=> []
push!(rate_constant_array,1.0);	# 105	addition_TGFB3_Ab: [] =([])=> TGFB3_Ab (reverse)
push!(rate_constant_array,0.1);	# 106	addition_ss_virus: ss_virus =([])=> []
push!(rate_constant_array,1.0);	# 107	addition_ss_virus: [] =([])=> ss_virus (reverse)
push!(rate_constant_array,0.1);	# 108	addition_smad_virus: smad_virus =([])=> []
push!(rate_constant_array,1.0);	# 109	addition_smad_virus: [] =([])=> smad_virus (reverse)
push!(rate_constant_array,0.1);	# 110	addition_lef1_virus: lef1_virus =([])=> []
push!(rate_constant_array,1.0);	# 111	addition_lef1_virus: [] =([])=> lef1_virus (reverse)
push!(rate_constant_array,1.0);	# 112	addition_TGFb12: [] =([])=> TGFB12
push!(rate_constant_array,0.1);	# 113	addition_TGFb12: TGFB12 =([])=> [] (reverse)
push!(rate_constant_array,1.0);	# 114	addition_VEGF: [] =([])=> VEGF
push!(rate_constant_array,0.1);	# 115	addition_VEGF: VEGF =([])=> [] (reverse)
push!(rate_constant_array,1.0);	# 116	addition_TGFB3: [] =([])=> TGFB3
push!(rate_constant_array,0.1);	# 117	addition_TGFB3: TGFB3 =([])=> [] (reverse)
push!(rate_constant_array,0.1);	# 118	degradation_mRNA_Smad: mRNA_Smad = []
push!(rate_constant_array,0.1);	# 119	degradation_mRNA_Ecadherin: mRNA_Ecadherin = []
push!(rate_constant_array,0.1);	# 120	degradation_mRNA_Bcatenin: mRNA_Bcatenin = []
push!(rate_constant_array,0.1);	# 121	degradation_mRNA_TCF4: mRNA_TCF4 = []
push!(rate_constant_array,0.1);	# 122	degradation_mRNA_TGFB3: mRNA_TGFB3 = []
push!(rate_constant_array,0.1);	# 123	degradation_mRNA_YREG1: mRNA_YREG1 = []
push!(rate_constant_array,0.1);	# 124	degradation_mRNA_LEF1: mRNA_LEF1 = []
push!(rate_constant_array,0.1);	# 125	degradation_mRNA_vimentin: mRNA_vimentin = []
push!(rate_constant_array,0.1);	# 126	degradation_mRNA_VEGFR: mRNA_VEGFR = []
push!(rate_constant_array,0.1);	# 127	degradation_mRNA_NFATc: mRNA_NFATc = []
push!(rate_constant_array,0.1);	# 128	degradation_mRNA_MAPK: mRNA_MAPK = []
push!(rate_constant_array,0.1);	# 129	degradation_mRNA_R: mRNA_R = []
push!(rate_constant_array,0.1);	# 130	degradation_mRNA_PLCy: mRNA_PLCy = []
push!(rate_constant_array,0.1);	# 131	degradation_mRNA_PKC: mRNA_PKC = []
push!(rate_constant_array,0.1);	# 132	degradation_mRNA_RAF: mRNA_RAF = []
push!(rate_constant_array,0.1);	# 133	degradation_mRNA_RAS: mRNA_RAS = []
push!(rate_constant_array,0.1);	# 134	degradation_mRNA_AP1_SP1: mRNA_AP1_SP1 = []
push!(rate_constant_array,0.1);	# 135	degradation_mRNA_Calcineurin: mRNA_Calcineurin = []
push!(rate_constant_array,0.1);	# 136	degradation_mRNA_SNAIL_SLUG: mRNA_SNAIL_SLUG = []
push!(rate_constant_array,0.1);	# 137	degradation_mRNA_tubulin: mRNA_tubulin = []
push!(rate_constant_array,0.1);	# 138	degradation_mRNA_GSK3_APC_AXIN: mRNA_GSK3_APC_AXIN = []
push!(rate_constant_array,0.1);	# 139	degradation_mRNA_PI3K: mRNA_PI3K = []
push!(rate_constant_array,1.0);	# 140	induction_gene_Smad: gene_Smad = mRNA_Smad
push!(rate_constant_array,1.0);	# 141	induction_gene_Ecadherin: gene_Ecadherin = mRNA_Ecadherin
push!(rate_constant_array,1.0);	# 142	induction_gene_Bcatenin: gene_Bcatenin = mRNA_Bcatenin
push!(rate_constant_array,1.0);	# 143	induction_gene_TCF4: gene_TCF4 = mRNA_TCF4
push!(rate_constant_array,1.0);	# 144	induction_gene_TGFB3: gene_TGFB3 = mRNA_TGFB3
push!(rate_constant_array,1.0);	# 145	induction_gene_YREG1: gene_YREG1 = mRNA_YREG1
push!(rate_constant_array,1.0);	# 146	induction_gene_LEF1: gene_LEF1 = mRNA_LEF1
push!(rate_constant_array,1.0);	# 147	induction_gene_R: gene_R = mRNA_R
push!(rate_constant_array,1.0);	# 148	induction_gene_vimentin: gene_vimentin = mRNA_vimentin
push!(rate_constant_array,1.0);	# 149	induction_gene_VEGFR: gene_VEGFR = mRNA_VEGFR
push!(rate_constant_array,1.0);	# 150	induction_gene_NFATc: gene_NFATc = mRNA_NFATc
push!(rate_constant_array,1.0);	# 151	induction_gene_PLCy: gene_PLCy = mRNA_PLCy
push!(rate_constant_array,1.0);	# 152	induction_gene_PKC: gene_PKC = mRNA_PKC
push!(rate_constant_array,1.0);	# 153	induction_gene_RAF: gene_RAF = mRNA_RAF
push!(rate_constant_array,1.0);	# 154	induction_gene_RAS: gene_RAS = mRNA_RAS
push!(rate_constant_array,1.0);	# 155	induction_gene_MAPK: gene_MAPK = mRNA_MAPK
push!(rate_constant_array,1.0);	# 156	induction_gene_AP1_SP1: gene_AP1_SP1 = mRNA_AP1_SP1
push!(rate_constant_array,1.0);	# 157	induction_gene_tubulin: gene_tubulin = mRNA_tubulin
push!(rate_constant_array,1.0);	# 158	induction_gene_SNAIL_SLUG: gene_SNAIL_SLUG = mRNA_SNAIL_SLUG
push!(rate_constant_array,1.0);	# 159	induction_gene_Calcineurin: gene_Calcineurin = mRNA_Calcineurin
push!(rate_constant_array,1.0);	# 160	induction_gene_GSK3_APC_AXIN: gene_GSK3_APC_AXIN = mRNA_GSK3_APC_AXIN
push!(rate_constant_array,1.0);	# 161	induction_gene_PI3K: gene_PI3K = mRNA_PI3K
push!(rate_constant_array,1.0);	# 162	translation_mRNA_Smad: mRNA_Smad = Smad
push!(rate_constant_array,1.0);	# 163	translation_mRNA_Ecadherin: mRNA_Ecadherin = Ecadherin
push!(rate_constant_array,1.0);	# 164	translation_mRNA_Bcatenin: mRNA_Bcatenin = Bcatenin
push!(rate_constant_array,1.0);	# 165	translation_mRNA_TCF4: mRNA_TCF4 = TCF4
push!(rate_constant_array,1.0);	# 166	translation_mRNA_TGFB3: mRNA_TGFB3 = TGFB3
push!(rate_constant_array,1.0);	# 167	translation_mRNA_YREG1: mRNA_YREG1 = YREG1
push!(rate_constant_array,1.0);	# 168	translation_mRNA_LEF1: mRNA_LEF1 = LEF1
push!(rate_constant_array,1.0);	# 169	translation_mRNA_vimentin: mRNA_vimentin = vimentin
push!(rate_constant_array,1.0);	# 170	translation_mRNA_VEGFR: mRNA_VEGFR = VEGFR
push!(rate_constant_array,1.0);	# 171	translation_mRNA_NFATc: mRNA_NFATc = NFATc
push!(rate_constant_array,1.0);	# 172	translation_mRNA_MAPK: mRNA_MAPK = MAPK
push!(rate_constant_array,1.0);	# 173	translation_mRNA_R: mRNA_R = R
push!(rate_constant_array,1.0);	# 174	translation_mRNA_PLCy: mRNA_PLCy = PLCy
push!(rate_constant_array,1.0);	# 175	translation_mRNA_PKC: mRNA_PKC = PKC
push!(rate_constant_array,1.0);	# 176	translation_mRNA_RAF: mRNA_RAF = RAF
push!(rate_constant_array,1.0);	# 177	translation_mRNA_RAS: mRNA_RAS = RAS
push!(rate_constant_array,1.0);	# 178	translation_mRNA_AP1_SP1: mRNA_AP1_SP1 = AP1_SP1
push!(rate_constant_array,1.0);	# 179	translation_mRNA_Calcineurin: mRNA_Calcineurin = Calcineurin
push!(rate_constant_array,1.0);	# 180	translation_mRNA_SNAIL_SLUG: mRNA_SNAIL_SLUG = SNAIL_SLUG
push!(rate_constant_array,1.0);	# 181	translation_mRNA_tubulin: mRNA_tubulin = tubulin
push!(rate_constant_array,1.0);	# 182	translation_mRNA_GSK3_APC_AXIN: mRNA_GSK3_APC_AXIN = GSK3_APC_AXIN
push!(rate_constant_array,1.0);	# 183	translation_mRNA_PI3K: mRNA_PI3K = PI3K

# Formulate the saturation constant array -
saturation_constant_array = zeros(NREACTIONS,NSPECIES);
saturation_constant_array[1,88] = 1.0;	# TGFB12_binding: TGFB12+R =([])=> TGFB12_R	 species: TGFB12
saturation_constant_array[1,87] = 1.0;	# TGFB12_binding: TGFB12+R =([])=> TGFB12_R	 species: R
saturation_constant_array[2,91] = 1.0;	# TGFB12_binding: TGFB12_R =([])=> TGFB12+R (reverse)	 species: TGFB12_R
saturation_constant_array[3,89] = 1.0;	# TGFB3_binding: TGFB3+R =([])=> TGFB3_R	 species: TGFB3
saturation_constant_array[3,87] = 1.0;	# TGFB3_binding: TGFB3+R =([])=> TGFB3_R	 species: R
saturation_constant_array[4,92] = 1.0;	# TGFB3_binding: TGFB3_R =([])=> TGFB3+R (reverse)	 species: TGFB3_R
saturation_constant_array[5,90] = 1.0;	# VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR	 species: VEGF
saturation_constant_array[5,93] = 1.0;	# VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR	 species: VEGFR
saturation_constant_array[6,94] = 1.0;	# VEGF_binding: VEGF_VEGFR =([])=> VEGF+VEGFR (reverse)	 species: VEGF_VEGFR
saturation_constant_array[7,27] = 1.0;	# VEGF_VEGFR_PLCy: PLCy =(VEGF_VEGFR)=> Active_PLCy	 species: PLCy
saturation_constant_array[8,29] = 1.0;	# PLCy_activation_Calcineurin: Calcineurin =(Active_PLCy)=> Active_Calcineurin	 species: Calcineurin
saturation_constant_array[9,31] = 1.0;	# Calcineurin_activation_NFATc: NFATc =(Active_Calcineurin)=> Active_NFATc	 species: NFATc
saturation_constant_array[10,33] = 1.0;	# PLCy_activatation_PKC: PKC =(Active_PLCy)=> Active_PKC	 species: PKC
saturation_constant_array[11,35] = 1.0;	# PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P	 species: RAF
saturation_constant_array[12,37] = 1.0;	# TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad	 species: Smad
saturation_constant_array[13,37] = 1.0;	# TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad	 species: Smad
saturation_constant_array[14,39] = 1.0;	# TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P	 species: PI3K
saturation_constant_array[15,41] = 1.0;	# PI3K_phosphorylation_GSK3: GSK3_APC_AXIN =(PI3K_P)=> GSK3_APC_AXIN_P	 species: GSK3_APC_AXIN
saturation_constant_array[16,43] = 1.0;	# TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP	 species: RAS
saturation_constant_array[17,35] = 1.0;	# RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P	 species: RAF
saturation_constant_array[18,45] = 1.0;	# RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK	 species: MAPK
saturation_constant_array[19,37] = 1.0;	# MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad	 species: Smad
saturation_constant_array[20,47] = 1.0;	# MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P	 species: AP1_SP1
saturation_constant_array[21,49] = 2.0;	# Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin	 species: Ecadherin
saturation_constant_array[21,50] = 0.5;	# Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin	 species: Bcatenin
saturation_constant_array[22,51] = 1.0;	# Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)	 species: Ecadherin_Bcatenin
saturation_constant_array[23,41] = 2.0;	# Binding_GSK3_APC_AXIN: GSK3_APC_AXIN+Bcatenin =([])=> GSK3_APC_AXIN_Bcatenin	 species: GSK3_APC_AXIN
saturation_constant_array[23,50] = 0.5;	# Binding_GSK3_APC_AXIN: GSK3_APC_AXIN+Bcatenin =([])=> GSK3_APC_AXIN_Bcatenin	 species: Bcatenin
saturation_constant_array[24,52] = 1.0;	# Binding_GSK3_APC_AXIN: GSK3_APC_AXIN_Bcatenin =([])=> GSK3_APC_AXIN+Bcatenin (reverse)	 species: GSK3_APC_AXIN_Bcatenin
saturation_constant_array[25,53] = 0.5;	# TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4	 species: TCF4
saturation_constant_array[25,50] = 2.0;	# TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4	 species: Bcatenin
saturation_constant_array[26,54] = 1.0;	# TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)	 species: Active_TCF4
saturation_constant_array[27,55] = 1.0;	# LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1	 species: LEF1
saturation_constant_array[27,50] = 1.0;	# LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1	 species: Bcatenin
saturation_constant_array[28,56] = 1.0;	# LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)	 species: Active_LEF1
saturation_constant_array[29,38] = 1.0;	# deactivation_pSmad: pSmad =(PPase)=> Smad	 species: pSmad
saturation_constant_array[30,36] = 1.0;	# deactivation_RAF_P: RAF_P =(PPase)=> RAF	 species: RAF_P
saturation_constant_array[31,44] = 1.0;	# deactivation_RAS_GTP: RAS_GTP =([])=> RAS	 species: RAS_GTP
saturation_constant_array[32,46] = 1.0;	# deactivation_pMAPK: pMAPK =(PPase)=> MAPK	 species: pMAPK
saturation_constant_array[33,48] = 1.0;	# deactivation_AP1_SP1: AP1_SP1_P =(PPase)=> AP1_SP1	 species: AP1_SP1_P
saturation_constant_array[34,28] = 1.0;	# deactivation_Active_PLCy: Active_PLCy =([])=> PLCy	 species: Active_PLCy
saturation_constant_array[35,34] = 1.0;	# deactivation_Active_PKC: Active_PKC =([])=> PKC	 species: Active_PKC
saturation_constant_array[36,30] = 1.0;	# deactivation_Active_Calcineurin: Active_Calcineurin =([])=> Calcineurin	 species: Active_Calcineurin
saturation_constant_array[37,32] = 1.0;	# deactivation_Active_NFATc: Active_NFATc =([])=> NFATc	 species: Active_NFATc
saturation_constant_array[38,42] = 1.0;	# deactivation_GSK3_APC_AXIN_P: GSK3_APC_AXIN_P =(PPase)=> GSK3_APC_AXIN	 species: GSK3_APC_AXIN_P

# Formulate control parameter array -
control_parameter_array = zeros(14,2);
# TGFB3_inhibition: TGFB3_Ab inhibition; TGFB3_binding
control_parameter_array[1,1] = 0.1;	#	1 Gain: 	TGFB3_inhibition target: TGFB3_binding actor: TGFB3_Ab type: inhibition
control_parameter_array[1,2] = 1.0;	#	1 Order: 	TGFB3_inhibition target: TGFB3_binding actor: TGFB3_Ab type: inhibition

# activation_Ecadherin_by_Active_NFATc: Active_NFATc induction induction_gene_Ecadherin
control_parameter_array[2,1] = 0.1;	#	2 Gain: 	activation_Ecadherin_by_Active_NFATc target: induction_gene_Ecadherin actor: Active_NFATc type: induction
control_parameter_array[2,2] = 1.0;	#	2 Order: 	activation_Ecadherin_by_Active_NFATc target: induction_gene_Ecadherin actor: Active_NFATc type: induction

# repression_Ecadherin_by_pSmad: pSmad repression induction_gene_Ecadherin
control_parameter_array[3,1] = 0.1;	#	3 Gain: 	repression_Ecadherin_by_pSmad target: induction_gene_Ecadherin actor: pSmad type: repression
control_parameter_array[3,2] = 1.0;	#	3 Order: 	repression_Ecadherin_by_pSmad target: induction_gene_Ecadherin actor: pSmad type: repression

# repression_Ecadherin_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_Ecadherin
control_parameter_array[4,1] = 0.1;	#	4 Gain: 	repression_Ecadherin_by_SNAIL_SLUG target: induction_gene_Ecadherin actor: SNAIL_SLUG type: repression
control_parameter_array[4,2] = 1.0;	#	4 Order: 	repression_Ecadherin_by_SNAIL_SLUG target: induction_gene_Ecadherin actor: SNAIL_SLUG type: repression

# repression_Ecadherin_by_Active_LEF1: Active_LEF1 repression induction_gene_Ecadherin
control_parameter_array[5,1] = 0.1;	#	5 Gain: 	repression_Ecadherin_by_Active_LEF1 target: induction_gene_Ecadherin actor: Active_LEF1 type: repression
control_parameter_array[5,2] = 1.0;	#	5 Order: 	repression_Ecadherin_by_Active_LEF1 target: induction_gene_Ecadherin actor: Active_LEF1 type: repression

# activation_TGFB3_by_TCF4: Active_TCF4 induction induction_gene_TGFB3
control_parameter_array[6,1] = 0.1;	#	6 Gain: 	activation_TGFB3_by_TCF4 target: induction_gene_TGFB3 actor: Active_TCF4 type: induction
control_parameter_array[6,2] = 1.0;	#	6 Order: 	activation_TGFB3_by_TCF4 target: induction_gene_TGFB3 actor: Active_TCF4 type: induction

# repression_YREG1_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_YREG1
control_parameter_array[7,1] = 0.1;	#	7 Gain: 	repression_YREG1_by_SNAIL_SLUG target: induction_gene_YREG1 actor: SNAIL_SLUG type: repression
control_parameter_array[7,2] = 1.0;	#	7 Order: 	repression_YREG1_by_SNAIL_SLUG target: induction_gene_YREG1 actor: SNAIL_SLUG type: repression

# activation_LEF1_by_pSmad: pSmad induction induction_gene_LEF1
control_parameter_array[8,1] = 0.1;	#	8 Gain: 	activation_LEF1_by_pSmad target: induction_gene_LEF1 actor: pSmad type: induction
control_parameter_array[8,2] = 1.0;	#	8 Order: 	activation_LEF1_by_pSmad target: induction_gene_LEF1 actor: pSmad type: induction

# repression_LEF1_by_YREG1: YREG1 repression induction_gene_LEF1
control_parameter_array[9,1] = 0.1;	#	9 Gain: 	repression_LEF1_by_YREG1 target: induction_gene_LEF1 actor: YREG1 type: repression
control_parameter_array[9,2] = 1.0;	#	9 Order: 	repression_LEF1_by_YREG1 target: induction_gene_LEF1 actor: YREG1 type: repression

# activation_vimentin_by_Active_LEF1: Active_LEF1 induction induction_gene_vimentin
control_parameter_array[10,1] = 0.1;	#	10 Gain: 	activation_vimentin_by_Active_LEF1 target: induction_gene_vimentin actor: Active_LEF1 type: induction
control_parameter_array[10,2] = 1.0;	#	10 Order: 	activation_vimentin_by_Active_LEF1 target: induction_gene_vimentin actor: Active_LEF1 type: induction

# activation_vimentin_by_AP1_SP1_P: AP1_SP1_P induction induction_gene_vimentin
control_parameter_array[11,1] = 0.1;	#	11 Gain: 	activation_vimentin_by_AP1_SP1_P target: induction_gene_vimentin actor: AP1_SP1_P type: induction
control_parameter_array[11,2] = 1.0;	#	11 Order: 	activation_vimentin_by_AP1_SP1_P target: induction_gene_vimentin actor: AP1_SP1_P type: induction

# activation_SNAIL_SLUG_by_AP1_SP1_P: AP1_SP1_P induction induction_gene_SNAIL_SLUG
control_parameter_array[12,1] = 0.1;	#	12 Gain: 	activation_SNAIL_SLUG_by_AP1_SP1_P target: induction_gene_SNAIL_SLUG actor: AP1_SP1_P type: induction
control_parameter_array[12,2] = 1.0;	#	12 Order: 	activation_SNAIL_SLUG_by_AP1_SP1_P target: induction_gene_SNAIL_SLUG actor: AP1_SP1_P type: induction

# activation_SNAIL_SLUG_by_ss_virus: ss_virus induction induction_gene_SNAIL_SLUG
control_parameter_array[13,1] = 0.1;	#	13 Gain: 	activation_SNAIL_SLUG_by_ss_virus target: induction_gene_SNAIL_SLUG actor: ss_virus type: induction
control_parameter_array[13,2] = 1.0;	#	13 Order: 	activation_SNAIL_SLUG_by_ss_virus target: induction_gene_SNAIL_SLUG actor: ss_virus type: induction

# Dilution selection array -
dilution_selection_array = ones(Int,NSPECIES)
dilution_selection_array[1] = 0.0;	#	gene_Smad
dilution_selection_array[2] = 0.0;	#	gene_Ecadherin
dilution_selection_array[3] = 0.0;	#	gene_Bcatenin
dilution_selection_array[4] = 0.0;	#	gene_PLCy
dilution_selection_array[5] = 0.0;	#	gene_PKC
dilution_selection_array[6] = 0.0;	#	gene_RAF
dilution_selection_array[7] = 0.0;	#	gene_RAS
dilution_selection_array[8] = 0.0;	#	gene_AP1_SP1
dilution_selection_array[9] = 0.0;	#	gene_tubulin
dilution_selection_array[10] = 0.0;	#	gene_SNAIL_SLUG
dilution_selection_array[11] = 0.0;	#	gene_Calcineurin
dilution_selection_array[12] = 0.0;	#	gene_GSK3_APC_AXIN
dilution_selection_array[13] = 0.0;	#	gene_PI3K
dilution_selection_array[69] = 0.0;	#	gene_TCF4
dilution_selection_array[70] = 0.0;	#	gene_TGFB3
dilution_selection_array[71] = 0.0;	#	gene_YREG1
dilution_selection_array[72] = 0.0;	#	gene_LEF1
dilution_selection_array[73] = 0.0;	#	gene_vimentin
dilution_selection_array[74] = 0.0;	#	gene_R
dilution_selection_array[75] = 0.0;	#	gene_VEGFR
dilution_selection_array[76] = 0.0;	#	gene_NFATc
dilution_selection_array[77] = 0.0;	#	gene_MAPK

# Set the maximum specific growth rate -
maximum_specific_growth_rate = 0.5;

# ---------------------------- DO NOT EDIT BELOW THIS LINE -------------------------- #
data_dictionary = Dict();
data_dictionary["STOICHIOMETRIC_MATRIX"] = S;
data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;
data_dictionary["SATURATION_CONSTANT_ARRAY"] = saturation_constant_array;
data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
data_dictionary["TIME_CONSTANT_ARRAY"] = time_constant_array;
data_dictionary["CONTROL_PARAMETER_ARRAY"] = control_parameter_array;
data_dictionary["MAXIMUM_SPECIFIC_GROWTH_RATE"] = maximum_specific_growth_rate;
data_dictionary["DILUTION_SELECTION_ARRAY"] = dilution_selection_array;
# ----------------------------------------------------------------------------------- #
return data_dictionary;
end