From 712da27cbd608ec9c0848ab6ae50f04d87f9f7a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Peter=20Br=C3=A4uer?= Date: Fri, 26 Oct 2018 12:10:11 +0100 Subject: [PATCH 1/4] Julia 1.0 updates Updates in function fit_j in fitTUV.jl to latest LsqFit version replacing estimate_errors with margin_error. Adding missing broadcast dot calculation of ss_tot. Additionally return fit data and errors from function j_oldpars in MCMphotolysis. Comment out the main script in the package. Remove Base packages from Project.toml. --- Project.toml | 3 --- src/MCMphotolysis.jl | 12 ++++++++---- src/fitTUV.jl | 5 ++--- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 5a6ca3a..ca3af3b 100644 --- a/Project.toml +++ b/Project.toml @@ -8,15 +8,12 @@ Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Juno = "e5e0dc1b-0480-54bc-9374-aad01c23163d" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" OhMyREPL = "5fb14364-9ced-5910-84b2-373655c76a03" OptimBase = "87e2bd06-a317-5318-96d9-3ecbac512eee" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" filehandling = "4b8ed5e0-8069-11e8-094b-3f5425e6acdd" diff --git a/src/MCMphotolysis.jl b/src/MCMphotolysis.jl index 7f684bb..bf9efda 100644 --- a/src/MCMphotolysis.jl +++ b/src/MCMphotolysis.jl @@ -55,17 +55,21 @@ function j_oldpars(scen) plot_jold(sza,χ,jvals,magnitude,fit,time,iofolder,scen) wrt_params(names(jvals),fit,magnitude,sigma,rmse,R2,iofolder,time) - # Return solar zenith angles, j values and fit data - return (sza, χ), jvals, fit + # Write data to output text file + =# - return (sza, χ), jvals + return (sza, χ), jvals, fit, (sigma, rmse, R2) end #function j_parameters end # module MCMphotolysis +#= +systime = Dates.now() -systime = now() println("load data...") scen = "testM4" ifile, iofolder = setup_files(scen) jvals, sza, χ = readTUV(ifile) + +scen = Juno.input("what? ") +=# diff --git a/src/fitTUV.jl b/src/fitTUV.jl index 2b6a503..559a24b 100644 --- a/src/fitTUV.jl +++ b/src/fitTUV.jl @@ -22,10 +22,10 @@ function fit_jold(jvals,sza) # Derive fit push!(fit, curve_fit(jold, sza, jvals[i], p0)) # Derive sigma with 95% confidence - push!(sigma, estimate_errors(fit[i],0.95)) + push!(sigma, margin_error(fit[i],0.05)) # Calculate statistical data for RMSE and R^2 ss_err = sum(fit[i].resid.^2) - ss_tot = sum((jvals[i]-mean(jvals[i])).^2) + ss_tot = sum((jvals[i].-mean(jvals[i])).^2) # RMSE push!(rmse, √(ss_err/fit[i].dof)) # R^2 @@ -33,7 +33,6 @@ function fit_jold(jvals,sza) end return fit, sigma, rmse, R2 - end #function fit_j From e0f9f3a2004e152fc57593806fdbfdf477814979 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Peter=20Br=C3=A4uer?= Date: Mon, 12 Nov 2018 14:55:02 +0000 Subject: [PATCH 2/4] Fix plotting with PyPlot to plot parameterisations Cleaning up and updating packages in the Project.toml and Manifest.toml. Initial README with installation instructions and some trouble shooting. Separating script and module code. The MCMphotolysis module is defined in MCMphotolysis.jl, script to call the Module found in a new file MCMparameters.jl. In MCMphotolysis.jl, updating package calls, fixing the PyCall of the PdfPages backend by importing it with the modules init functions (thanks to @kilicomu). New flag `output` for function j_oldpars to plot parameterisations and write parameters to file. New include file output.jl, to use PyPlot to plot parameterisations. In rdfiles.jl, use test_file from filehandling to read TUV files, which asks for user input, if a file name doesn't exist rather than stopping the script. Only create output folder and ask for overwrite, if output flag is set to true. --- Manifest.toml | 82 ++++++++++++++++++++++++------------------ Project.toml | 11 +++--- README.md | 44 +++++++++++++++++++++++ src/MCMparameters.jl | 6 ++++ src/MCMphotolysis.jl | 62 +++++++++++++++----------------- src/output.jl | 85 ++++++++++++++++++++++++++++++++++++++++++++ src/rdfiles.jl | 10 +++--- 7 files changed, 221 insertions(+), 79 deletions(-) create mode 100644 src/MCMparameters.jl create mode 100644 src/output.jl diff --git a/Manifest.toml b/Manifest.toml index 3e1cac8..0636e4e 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -32,10 +32,10 @@ uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597" version = "0.4.0" [[CodecZlib]] -deps = ["BinaryProvider", "Libdl", "Pkg", "Test", "TranscodingStreams"] -git-tree-sha1 = "83cb3d65c37ea1364c2d5bf7bcea41843ba645dc" +deps = ["BinaryProvider", "Libdl", "Test", "TranscodingStreams"] +git-tree-sha1 = "e3df104c84dfc108f0ca203fd7f5bbdc98641ae9" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.5.0" +version = "0.5.1" [[ColorTypes]] deps = ["FixedPointNumbers", "Random", "Test"] @@ -44,7 +44,7 @@ uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" version = "0.7.5" [[Colors]] -deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Pkg", "Printf", "Reexport", "Test"] +deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Printf", "Reexport", "Test"] git-tree-sha1 = "9f0a0210450acb91c730b730a994f8eef1d3d543" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" version = "0.9.5" @@ -68,25 +68,25 @@ uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" version = "1.1.1" [[Crayons]] -deps = ["Pkg", "Test"] +deps = ["Test"] git-tree-sha1 = "3017c662a988bcb8a3f43306a793617c6524d476" uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "1.0.0" [[DataFrames]] -deps = ["CategoricalArrays", "CodecZlib", "Compat", "DataStreams", "Dates", "InteractiveUtils", "IteratorInterfaceExtensions", "LinearAlgebra", "Missings", "Pkg", "Printf", "Random", "Reexport", "SortingAlgorithms", "Statistics", "StatsBase", "TableTraits", "Tables", "Test", "TranscodingStreams", "Unicode", "WeakRefStrings"] +deps = ["CategoricalArrays", "CodecZlib", "Compat", "DataStreams", "Dates", "InteractiveUtils", "IteratorInterfaceExtensions", "LinearAlgebra", "Missings", "Printf", "Random", "Reexport", "SortingAlgorithms", "Statistics", "StatsBase", "TableTraits", "Tables", "Test", "TranscodingStreams", "Unicode", "WeakRefStrings"] git-tree-sha1 = "ad34fefb72b18a8dd5c17fab9089d11111b61935" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" version = "0.14.1" [[DataStreams]] -deps = ["Dates", "Missings", "Pkg", "Test", "WeakRefStrings"] +deps = ["Dates", "Missings", "Test", "WeakRefStrings"] git-tree-sha1 = "69c72a1beb4fc79490c361635664e13c8e4a9548" uuid = "9a8bc11e-79be-5b39-94d7-1ccc349a1a85" version = "0.4.1" [[DataStructures]] -deps = ["InteractiveUtils", "OrderedCollections", "REPL", "Random", "Serialization", "Test"] +deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"] git-tree-sha1 = "8fc6e166e24fda04b2b648d4260cdad241788c54" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" version = "0.14.0" @@ -101,9 +101,9 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" [[DiffEqDiffTools]] deps = ["LinearAlgebra", "Test"] -git-tree-sha1 = "5d220295e07640ae3c096c3511fcc15081e61f80" +git-tree-sha1 = "67700c9fc82033ec68a145bc650f6b9debdf9103" uuid = "01453d9d-ee7c-5054-8395-0335cb756afa" -version = "0.7.0" +version = "0.7.1" [[DiffResults]] deps = ["Compat", "StaticArrays"] @@ -122,7 +122,7 @@ deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[Distributions]] -deps = ["Distributed", "LinearAlgebra", "PDMats", "Pkg", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"] +deps = ["Distributed", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"] git-tree-sha1 = "c24e9b6500c037673f0241a2783472b8c3d080c7" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" version = "0.16.4" @@ -131,16 +131,16 @@ version = "0.16.4" uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" [[FixedPointNumbers]] -deps = ["Pkg", "Test"] +deps = ["Test"] git-tree-sha1 = "b8045033701c3b10bf2324d7203404be7aef88ba" uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" version = "0.5.3" [[ForwardDiff]] -deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Pkg", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"] -git-tree-sha1 = "d8f3e0f19d0d546aa92eb1cd67cd3e515768d9f7" +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"] +git-tree-sha1 = "b91250044374764e7c29af59a774c4b8d6100b6e" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.0" +version = "0.10.1" [[Future]] deps = ["Random"] @@ -157,7 +157,7 @@ uuid = "82899510-4779-5014-852e-03e436cf321d" version = "0.1.1" [[JSON]] -deps = ["Dates", "Distributed", "Mmap", "Pkg", "Sockets", "Test", "Unicode"] +deps = ["Dates", "Distributed", "Mmap", "Sockets", "Test", "Unicode"] git-tree-sha1 = "fec8e4d433072731466d37ed0061b3ba7f70eeb9" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.19.0" @@ -188,7 +188,7 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[LsqFit]] -deps = ["Calculus", "Distributions", "LinearAlgebra", "OptimBase", "Pkg", "Random", "Test"] +deps = ["Calculus", "Distributions", "LinearAlgebra", "OptimBase", "Random", "Test"] git-tree-sha1 = "5c1b5ab85e120c6d1734e7f5670aac22fe69270c" uuid = "2fda8390-95c7-5789-9bda-21331edee243" version = "0.6.0" @@ -205,9 +205,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[Media]] deps = ["MacroTools", "Test"] -git-tree-sha1 = "9f390271c9a43dcbe908a10b5b9632cf58cbab5b" +git-tree-sha1 = "75a54abd10709c01f1b86b84ec225d26e840ed58" uuid = "e89f7d12-3494-54d1-8411-f7d8b9ae1f27" -version = "0.4.1" +version = "0.5.0" [[Missings]] deps = ["Dates", "InteractiveUtils", "SparseArrays", "Test"] @@ -219,7 +219,7 @@ version = "0.3.1" uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[NLSolversBase]] -deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff", "LinearAlgebra", "Pkg", "Random", "SparseArrays", "Test"] +deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff", "LinearAlgebra", "Random", "SparseArrays", "Test"] git-tree-sha1 = "ebfb2e96970151753575b9c4d31d47e5ae8382a5" uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" version = "7.1.1" @@ -231,19 +231,19 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "0.3.2" [[OhMyREPL]] -deps = ["Crayons", "InteractiveUtils", "Markdown", "Pkg", "Printf", "REPL", "Test", "Tokenize"] +deps = ["Crayons", "Markdown", "Pkg", "Printf", "REPL", "Test", "Tokenize"] git-tree-sha1 = "e00d5394d110afe279101ffe10cebd11eaedcb8a" uuid = "5fb14364-9ced-5910-84b2-373655c76a03" version = "0.3.0" [[OptimBase]] -deps = ["Compat", "NLSolversBase", "Pkg", "Printf", "Reexport", "Test"] +deps = ["Compat", "NLSolversBase", "Printf", "Reexport", "Test"] git-tree-sha1 = "92667ab46a66ad502ec3044f65c41ea68b2e0e9c" uuid = "87e2bd06-a317-5318-96d9-3ecbac512eee" version = "2.0.0" [[OrderedCollections]] -deps = ["Pkg", "Random", "Serialization", "Test"] +deps = ["Random", "Serialization", "Test"] git-tree-sha1 = "85619a3f3e17bb4761fe1b1fd47f0e979f964d5b" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.0.2" @@ -254,6 +254,12 @@ git-tree-sha1 = "9e3e7a5c9b8cfdba8c01a1bcae38fe0144936fda" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" version = "0.9.5" +[[Parameters]] +deps = ["Markdown", "OrderedCollections", "REPL", "Test"] +git-tree-sha1 = "40f540ec96e50c0b2b9efdb11b5e4d0c63f90923" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.10.1" + [[Pkg]] deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" @@ -317,7 +323,7 @@ uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" version = "0.7.12" [[Rmath]] -deps = ["BinaryProvider", "Libdl", "Pkg", "Random", "Statistics", "Test"] +deps = ["BinaryProvider", "Libdl", "Random", "Statistics", "Test"] git-tree-sha1 = "9a6c758cdf73036c3239b0afbea790def1dabff9" uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" version = "0.5.0" @@ -347,12 +353,12 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SpecialFunctions]] deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"] -git-tree-sha1 = "c35c9c76008babf4d658060fc64aeb369a41e7bd" +git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "0.7.1" +version = "0.7.2" [[StaticArrays]] -deps = ["InteractiveUtils", "LinearAlgebra", "Pkg", "Random", "Statistics", "Test"] +deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"] git-tree-sha1 = "ebc5c2a27d91d5ec611a9861168182e2168effd3" uuid = "90137ffa-7385-5640-81b9-e52037218182" version = "0.9.2" @@ -368,7 +374,7 @@ uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" version = "0.25.0" [[StatsFuns]] -deps = ["Pkg", "Rmath", "SpecialFunctions", "Test"] +deps = ["Rmath", "SpecialFunctions", "Test"] git-tree-sha1 = "d14bb7b03defd2deaa5675646f6783089e0556f0" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" version = "0.7.0" @@ -384,10 +390,10 @@ uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" version = "0.4.0" [[Tables]] -deps = ["Pkg", "Requires", "Test"] -git-tree-sha1 = "be3a3a28a1e6c4fdd838bdf21e796f7a592c565c" +deps = ["Requires", "Test"] +git-tree-sha1 = "c7fb447deab835fa70ce6717e78c68b0f466a42c" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "0.1.10" +version = "0.1.11" [[Test]] deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] @@ -400,7 +406,7 @@ uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624" version = "0.5.2" [[TranscodingStreams]] -deps = ["DelimitedFiles", "Pkg", "Random", "Test"] +deps = ["Pkg", "Random", "Test"] git-tree-sha1 = "a34a2d588e2d2825602bf14a24216d5c8b0921ec" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" version = "0.8.1" @@ -431,9 +437,17 @@ uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5" version = "0.5.3" [[filehandling]] -deps = ["DataFrames", "Revise"] -git-tree-sha1 = "9519bbba686f41600e696ccd06b29bba3235f02a" +deps = ["DataFrames", "Juno", "Revise"] +git-tree-sha1 = "1b9ed9e6a80d07ebdd8b903c85f264f6a374b8e5" repo-rev = "master" repo-url = "https://github.com/pb866/filehandling.git" uuid = "4b8ed5e0-8069-11e8-094b-3f5425e6acdd" +version = "0.1.2" + +[[pyp]] +deps = ["DataFrames", "Dates", "LaTeXStrings", "LinearAlgebra", "OhMyREPL", "Parameters", "PyPlot", "filehandling"] +git-tree-sha1 = "85d1e7c993f1b3e452f2fde19780f78ea1f07310" +repo-rev = "master" +repo-url = "https://github.com/pb866/pyp.git" +uuid = "1f59166c-d923-11e8-01b9-e713fdda0351" version = "0.1.0" diff --git a/Project.toml b/Project.toml index ca3af3b..6d3b239 100644 --- a/Project.toml +++ b/Project.toml @@ -4,16 +4,15 @@ authors = ["Peter Bräuer "] version = "0.1.0" [deps] -Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" -Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" -Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Juno = "e5e0dc1b-0480-54bc-9374-aad01c23163d" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" OhMyREPL = "5fb14364-9ced-5910-84b2-373655c76a03" -OptimBase = "87e2bd06-a317-5318-96d9-3ecbac512eee" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" +PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" filehandling = "4b8ed5e0-8069-11e8-094b-3f5425e6acdd" +pyp = "1f59166c-d923-11e8-01b9-e713fdda0351" diff --git a/README.md b/README.md index b8ed2dd..9c34a7f 100644 --- a/README.md +++ b/README.md @@ -1 +1,45 @@ MCMphotolysis +============= + +A Julia package to retrieve updated MCM photolysis parameterisations including +a dependence on the overlying ozone column from TUV output files. + + +Installation +------------ + +Install the Package by adding it to your environment using `Pkg` + +```julia +using Pkg +Pkg.add("https://github.com/pb866/MCMphotolysis.git") +Pkg.precomile() +``` + +or go to the package manager typing `]` in the julia prompt and in the +package manager the following commands: + +```julia +add https://github.com/pb866/MCMphotolysis.git +precompile +``` + +### Trouble shooting of common errors + +If you get an error message, follow the the instructions of the error message, e.g. + +```julia +Pkg.build("CodecZlib") +``` + +If PyPlot crashes, try running Julia with the system python rather than +the miniconda python version by rebuilding python with: + +```julia +using Pkg +ENV["PYTHON"] = "path/to/python" +Pkg.build("PyPlot") +``` + +You can get the system python version by typing `which python` in the terminal +and copying the output. diff --git a/src/MCMparameters.jl b/src/MCMparameters.jl new file mode 100644 index 0000000..1e993e6 --- /dev/null +++ b/src/MCMparameters.jl @@ -0,0 +1,6 @@ +using Pkg +Pkg.activate(".") +using MCMphotolysis + +jvals = j_oldpars("testM4", output = false) +jvals = j_oldpars("testM4") diff --git a/src/MCMphotolysis.jl b/src/MCMphotolysis.jl index bf9efda..d33a4d3 100644 --- a/src/MCMphotolysis.jl +++ b/src/MCMphotolysis.jl @@ -3,20 +3,22 @@ module MCMphotolysis # Source directory dir = Base.source_dir() -# Track changes during development -using Pkg -Pkg.activate(".") -using Revise - # Load Julia packages using Statistics +using LinearAlgebra using Juno: input #get terminal input inside or outside Atom -using Dates: now +using Dates using ProgressMeter using LsqFit +using PyPlot, PyCall +using DataFrames +# PyCall python imports +# @pyimport matplotlib.backends.backend_pdf as pdf +const pdf = PyNULL() # Load self-made packages using filehandling +import pyp # export public functions export j_oldpars @@ -24,6 +26,11 @@ export j_oldpars # Include outsourced functions include(joinpath(dir,"rdfiles.jl")) include(joinpath(dir,"fitTUV.jl")) +include(joinpath(dir,"output.jl")) + +function __init__() + copy!(pdf, pyimport("matplotlib.backends.backend_pdf")) +end """ j_oldpars(scen) @@ -35,41 +42,30 @@ to a pdf and print the j values and errors/statistical data to a text file in a Return solar zenith angle as tuple of arrays with deg/rad, the TUV _j_ value data, and the fit data (as array with entries for each reaction). """ -function j_oldpars(scen) +function j_oldpars(scen::String; output::Bool=true) # Initialise system time and output path/file name systime = now() # Read dataframe with j values from TUV output file println("load data...") - ifile, iofolder = setup_files(scen) - jvals, sza, χ = readTUV(ifile) + ifile, iofolder = setup_files(scen, output) + jvals = readTUV(ifile) - magnitude = [floor(log10(jvals[i][1])) for i = 1:length(jvals)] - for i = 1:length(magnitude) jvals[i] ./= 10^magnitude[i] end + magnitude = [floor(log10(jvals[:jvals][i][1])) for i = 1:length(jvals[:jvals])] + for i = 1:length(magnitude) jvals[:jvals][i] ./= 10^magnitude[i] end + jvals[:order] = magnitude # Derive parameterisations for j values - fit, sigma, rmse, R2 = fit_jold(jvals,χ) - - #= - # Generate output - plot_jold(sza,χ,jvals,magnitude,fit,time,iofolder,scen) - wrt_params(names(jvals),fit,magnitude,sigma,rmse,R2,iofolder,time) - - # Write data to output text file - - =# - return (sza, χ), jvals, fit, (sigma, rmse, R2) + fit, sigma, rmse, R2 = fit_jold(jvals[:jvals],jvals[:rad]) + jvals[:fit] = fit; jvals[:σ] = sigma, jvals[:RMSE] = rmse; jvals[:R2] = R2 + + # # Generate output + if output + # println("ϑ") + plot_jold(jvals,systime,iofolder,scen) + # wrt_params(names(jvals),fit,magnitude,sigma,rmse,R2,iofolder,time) + end + return jvals end #function j_parameters end # module MCMphotolysis - -#= -systime = Dates.now() - -println("load data...") -scen = "testM4" -ifile, iofolder = setup_files(scen) -jvals, sza, χ = readTUV(ifile) - -scen = Juno.input("what? ") -=# diff --git a/src/output.jl b/src/output.jl new file mode 100644 index 0000000..dfa75d7 --- /dev/null +++ b/src/output.jl @@ -0,0 +1,85 @@ +""" + plot_j(jvals, iofolder, fit, systime) + +Plot j values saved in dataframe jvals and parameterisations derived from parameters +in fit to file parameters.dat in the iofolder together with the time of creation (time). +""" +function plot_jold(jvals::Dict{Symbol,Any},systime::DateTime,iofolder::String,ifile::String) + + # Format titles with Latex using header of jvals + ptitle = beautify_chem(names(jvals[:jvals])) + + # Initialise array of plots and define x data + nj = length(jvals[:fit]) + # jplot = Array{PyCall.PyObject}(nj) + ofile = "$iofolder/$ifile.pdf" + opfile = pdf[:PdfPages](ofile) + + # Loop over all reactions + @showprogress 1 "plot data..." for i=1:nj + # define parameters + l = jvals[:fit][i].param[1] + m = jvals[:fit][i].param[2] + n = jvals[:fit][i].param[3] + # Calculate parameterised values + jpar = l.⋅cos.(jvals[:rad]).^m.⋅exp.(-n.⋅sec.(jvals[:rad])) + # Load TUV data and parameterisation for plotting + jplt = pyp.load_PlotData(DataFrame(χ=jvals[:deg], j=jvals[:jvals][i]), + label = "TUV data", pt = "s", lc = "black", lt = "None") + pplt = pyp.load_PlotData(DataFrame(x=jvals[:deg], y=jpar), lc ="red", lw=2, + label="MCM Parameterisaton") + # Plot TUV data + fig, ax = pyp.plot_data(jplt, pplt, ti = ptitle[i], + xlabel = "solar zenith angle χ", + ylabel = "j / 10\$^{$(Int(jvals[:order][i]))}\$ s\$^{-1}\$", + xlims=(0,90), ylims=(0,nothing), maj_xticks=15, min_xticks=5, ti_offset=-2, + legpos="lower left") + + # Plot time stamp + ax[:annotate]("created $(Dates.format(systime,"dd.mm.yyyy, HH:MM:SS"))", + xy=[0.02;0.22], xycoords="axes fraction", ha="left", va="bottom") + + # save plot to temporary png + opfile[:savefig](fig) + close() + end + + # compile pngs in single pdf and delete pngs + opfile[:close]() +end #function plot_j + + +""" + beautify_chem(reactions) + +Format `reactions` with LaTeX for nicer titles in plots. +""" +function beautify_chem(reactions::Vector{Symbol}) + + # Initialise output + chem = String[] + # Loop over reactions and reformat with LaTeX + for rxn in strip.(string.(reactions)) + println(rxn) + lstr = replace(rxn, "->" => "\\stackrel{h\\nu}{\\longrightarrow}") # format arrows + lstr = replace(lstr, " + hv" => "") # remove + hv (now above formatted arrows) + lstr = replace(lstr, r"([A-Za-z)])([0-9]+)" => s"\1_{\2}") # make numbers in chemical formulas subscripts: {\1}_{\2} + lstr = replace(lstr, r"\(([0-9]+)" => s"(^{\1}") # except for atom's energetic state + lstr = replace(lstr, "." => "^{.}") # raise radical dots + lstr = replace(lstr, r"([0-9]+) " => s"\1\\,") # insert halfspace between stoichiometric indices and species + lstr = replace(lstr, r"=" => "\\!=\\!") # no spaces between double bonds + lstr = replace(lstr, r"([A-Za-z0-9])-([A-Za-z])" => s"\1\\!-\\!\2") # format dashes + lstr = replace(lstr, r"C:" => "{\\ddot C}") # format biradical functions + lstr = replace(lstr, r"CH:" => "{\\ddot C}H") # format biradical functions + # define and correct exceptions: + lstr = replace(lstr, " " => "\\ ") # ensure spaces + # Ensure unslanted font + lstr = "\\mathrm{"*lstr*"}" + + # Save reformatted reaction to output + push!(chem, lstr) + end + + # Return reformatted output + return latexstring.(chem) +end #function beautify_chem diff --git a/src/rdfiles.jl b/src/rdfiles.jl index 146374c..e313548 100644 --- a/src/rdfiles.jl +++ b/src/rdfiles.jl @@ -4,7 +4,7 @@ From the scenario name (scen) used in TUV derive names for the TUV output file and the output folder. """ -function setup_files(scen) +function setup_files(scen, output) # Test script argument and ask for input file, if missing if isempty(scen) println("Enter name for TUV scenario") @@ -13,7 +13,7 @@ function setup_files(scen) # Create output folder iofolder = "./params_"*strip(scen) - try mkdir(iofolder) + if output try mkdir(iofolder) catch print("\033[95mFolder '$iofolder' already exists! Overwrite ") confirm = input("(\033[4mY\033[0m\033[95mes/\033[4mN\033[0m\033[95mo)?\033[0m ") @@ -21,13 +21,11 @@ function setup_files(scen) cd(iofolder); files = readdir(); [rm(file) for file in files ]; cd("..") else println("Programme aborted. Exit code '98'."); exit(98) end - end + end end # Define TUV file ifile = scen*".txt" - if !isfile(ifile) - println("$ifile does not exist. Script terminated. Exit code: 99."); exit(99) - end + ifile = test_file(ifile) # return file and folder names return ifile, iofolder From 9ebd2383d0fa5b431043dfb396a1fc83067b8d38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Peter=20Br=C3=A4uer?= Date: Mon, 12 Nov 2018 16:24:37 +0000 Subject: [PATCH 3/4] Write legacy parameters to output file, fix Package Printf added to Project.toml for formatted text file output. New function wrt_params to write legacy MCM parameters to text file in output.jl. Adding call of Printf and function wrt_params to MCMphotolysis.jl. FIX: Changing comma to semicolon before entry RMSE, when appending the jvals dictionary on l. 69 in MCMphotolysis.jl. --- Project.toml | 1 + src/MCMphotolysis.jl | 5 +++-- src/output.jl | 39 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 42 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 6d3b239..ebbcfd6 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ Juno = "e5e0dc1b-0480-54bc-9374-aad01c23163d" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" OhMyREPL = "5fb14364-9ced-5910-84b2-373655c76a03" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" diff --git a/src/MCMphotolysis.jl b/src/MCMphotolysis.jl index d33a4d3..27a6b3a 100644 --- a/src/MCMphotolysis.jl +++ b/src/MCMphotolysis.jl @@ -12,6 +12,7 @@ using ProgressMeter using LsqFit using PyPlot, PyCall using DataFrames +using Printf # PyCall python imports # @pyimport matplotlib.backends.backend_pdf as pdf const pdf = PyNULL() @@ -57,13 +58,13 @@ function j_oldpars(scen::String; output::Bool=true) # Derive parameterisations for j values fit, sigma, rmse, R2 = fit_jold(jvals[:jvals],jvals[:rad]) - jvals[:fit] = fit; jvals[:σ] = sigma, jvals[:RMSE] = rmse; jvals[:R2] = R2 + jvals[:fit] = fit; jvals[:σ] = sigma; jvals[:RMSE] = rmse; jvals[:R2] = R2 # # Generate output if output # println("ϑ") plot_jold(jvals,systime,iofolder,scen) - # wrt_params(names(jvals),fit,magnitude,sigma,rmse,R2,iofolder,time) + wrt_params(jvals,iofolder,systime) end return jvals end #function j_parameters diff --git a/src/output.jl b/src/output.jl index dfa75d7..498d4cf 100644 --- a/src/output.jl +++ b/src/output.jl @@ -11,7 +11,6 @@ function plot_jold(jvals::Dict{Symbol,Any},systime::DateTime,iofolder::String,if # Initialise array of plots and define x data nj = length(jvals[:fit]) - # jplot = Array{PyCall.PyObject}(nj) ofile = "$iofolder/$ifile.pdf" opfile = pdf[:PdfPages](ofile) @@ -49,6 +48,44 @@ function plot_jold(jvals::Dict{Symbol,Any},systime::DateTime,iofolder::String,if end #function plot_j +""" + wrt_params(rxn, fit, sigma, rmse, R2, iofolder, time) + +For each reaction (rxn), from param in fit and sigma, print parameters +and 95% confidence together with RMSE and R^2 to file 'parameters.dat' +in the designated output folder (iofolder). + +Print the time of creation (time) to the output file. +""" +function wrt_params(jvals, iofolder, systime) + + #transform dataframe column symbols to strings + rxn = string.(names(jvals[:jvals])) + + # Open output file + open("$iofolder/parameters.dat","w") do f + # Print header + println(f, + "Parameters and statistical data for parameterisation of photolysis processes") + println(f, + "in the Master Chemical Mechanism (MCM; http://mcm.york.ac.uk/) using:") + println(f, "\nj / s-1 = l·(cos(x))^m·exp(-n·sec(x))") + println(f, "\n\ncreated $(Dates.format(systime,"dd.mm.yyyy, HH:MM:SS"))") + println(f,"\n P a r a m e t e r s S t a t i s t i c s") + println(f," l / s-1 m n RMSE / s-1 R^2 Reaction") + + # Loop over reactions + for i = 1:length(jvals[:fit]) + # Print parameters, statistical data, and reaction label to output file + @printf(f,"(%6.3f±%.3f)e%d %.3f±%.3f %.3f±%.3f %.3e %.4f %s\n", + jvals[:fit][i].param[1], jvals[:σ][i][1], Int(jvals[:order][i]), + jvals[:fit][i].param[2], jvals[:σ][i][2], jvals[:fit][i].param[3], jvals[:σ][i][3], + jvals[:RMSE][i]⋅10^jvals[:order][i], jvals[:R2][i], rxn[i]) + end + end +end # function wrt_params + + """ beautify_chem(reactions) From 279a8315b7e16661110c3e913341cfc45e273cdc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Peter=20Br=C3=A4uer?= Date: Mon, 12 Nov 2018 17:08:00 +0000 Subject: [PATCH 4/4] Updated documentation Update README with usage of function j_oldpars. Add version history. Add Julia version ("0.7-") to Project.toml. Remove script MCMparameters.jl from src folder. Updated docstrings. --- Project.toml | 1 + README.md | 46 ++++++++++++++++++++++++++++++++++++++++++-- src/MCMparameters.jl | 6 ------ src/MCMphotolysis.jl | 20 +++++++++++++------ src/fitTUV.jl | 8 +++++--- src/output.jl | 17 ++++++++-------- src/rdfiles.jl | 7 +++++-- 7 files changed, 77 insertions(+), 28 deletions(-) delete mode 100644 src/MCMparameters.jl diff --git a/Project.toml b/Project.toml index ebbcfd6..dd8d572 100644 --- a/Project.toml +++ b/Project.toml @@ -2,6 +2,7 @@ name = "MCMphotolysis" uuid = "37f95c2c-7f9b-11e8-13ad-816be0d1e117" authors = ["Peter Bräuer "] version = "0.1.0" +julia = "0.7-" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" diff --git a/README.md b/README.md index 9c34a7f..e23081b 100644 --- a/README.md +++ b/README.md @@ -13,18 +13,49 @@ Install the Package by adding it to your environment using `Pkg` ```julia using Pkg Pkg.add("https://github.com/pb866/MCMphotolysis.git") +Pkg.instantiate() Pkg.precomile() ``` or go to the package manager typing `]` in the julia prompt and in the package manager the following commands: -```julia +``` add https://github.com/pb866/MCMphotolysis.git +instantiate precompile ``` -### Trouble shooting of common errors + +Usage +----- + +After installation run the following lines in the REPL or write a julia +script with the following lines and call it from the terminal. + +```julia +using MCMphotolysis +jvals = j_oldpars("") +``` + +If you haven't installed the package directly to your Julia default environment, +you need to activate the environment first. + +Call function `j_oldpars` handing over the scenario name, i.e. the name of +the TUV input file without the `.txt` file ending. If you only want the data +returned in the REPL without file output, set the optional keyword argument +`output = false`. +Function `j_oldpars` returns a dictionary with entries for a DataFrame with +the _j_ values (`:jvals`), where the `:order` of magnitude is stored in a +different array, arrays with the solor zenith angles (`:deg`/`:rad`), +an array with LsqFit output, statistical data (`:σ`, `:RMSE`, and `:R2`). +Moreover, a folder `params_` is created, where parameters are +printed to `parameters.dat` and parameterisations are visualised in +`.pdf`. + + +Trouble shooting of common errors +--------------------------------- If you get an error message, follow the the instructions of the error message, e.g. @@ -43,3 +74,14 @@ Pkg.build("PyPlot") You can get the system python version by typing `which python` in the terminal and copying the output. + + +Version history +=============== + +Version 0.1.0 +------------- + +- Function `j_oldpars` to write the legacy MCM parameterisations and statistical + data to a text file and plot the parameterisations in a pdf compiled in the + folder `params_` diff --git a/src/MCMparameters.jl b/src/MCMparameters.jl deleted file mode 100644 index 1e993e6..0000000 --- a/src/MCMparameters.jl +++ /dev/null @@ -1,6 +0,0 @@ -using Pkg -Pkg.activate(".") -using MCMphotolysis - -jvals = j_oldpars("testM4", output = false) -jvals = j_oldpars("testM4") diff --git a/src/MCMphotolysis.jl b/src/MCMphotolysis.jl index 27a6b3a..a53eab5 100644 --- a/src/MCMphotolysis.jl +++ b/src/MCMphotolysis.jl @@ -34,14 +34,22 @@ function __init__() end """ - j_oldpars(scen) + j_oldpars(scen::String; output::Bool=true) Read reactions from a TUV output file named `scen`.txt, derive MCM l,m,n-parameterisations -for photolysis, plot the _j_ values from TUV and the parameterisation for all reactions -to a pdf and print the j values and errors/statistical data to a text file in a folder -`params_`. -Return solar zenith angle as tuple of arrays with deg/rad, the TUV _j_ value data, -and the fit data (as array with entries for each reaction). +for photolysis and return a dictionary with the following entries: + +- `:jvals`: DataFrame with j values (devided by maximum order of magnitude, i.e. without the `e-...`) +- `:order`: Vector with order of magnitudes +- `:deg`/`:rad`: Vector with solar zenith angles in deg/rad +- `:fit`: Vector with LsqFit data +- `:σ`: Vector of Vectors with σ-values for the 95% confidence interval for the `l`, `m`, and `n` parameters +- `:RMSE`: Vector of root mean square errors +- `:R2`: Vector of correlation cofficients R² + +If output is set to `true` (_default_), _j_ values from TUV and the parameterisations +for all reactions are plotted to a pdf and _j_ values and errors/statistical data +are printed to a text file in a folder named `params_`. """ function j_oldpars(scen::String; output::Bool=true) # Initialise system time and output path/file name diff --git a/src/fitTUV.jl b/src/fitTUV.jl index 559a24b..2ff6884 100644 --- a/src/fitTUV.jl +++ b/src/fitTUV.jl @@ -3,10 +3,12 @@ ######################################### """ - fit_j(jvals) + fit_j(jvals, sza) -Derive MCM parameterisations for dataframe jvals with χ-dependent _j_ values. - j / s^-1 = l·cos^m(χ)·exp(-n·sec(χ)) +Derive MCM parameterisations for DataFrame `jvals` with `sza`-dependent _j_ values. + j / s^-1 = l·cos^m(sza)·exp(-n·sec(sza))] + +The Vector of solar zenith angles `sza` must be in rad. """ function fit_jold(jvals,sza) diff --git a/src/output.jl b/src/output.jl index 498d4cf..8d49502 100644 --- a/src/output.jl +++ b/src/output.jl @@ -1,8 +1,9 @@ """ - plot_j(jvals, iofolder, fit, systime) + plot_j(jvals::Dict{Symbol,Any},systime::DateTime,iofolder::String,ifile::String) -Plot j values saved in dataframe jvals and parameterisations derived from parameters -in fit to file parameters.dat in the iofolder together with the time of creation (time). +Plot j values saved in `jvals[:jvals]` and parameterisations derived from parameters +in `jvals[:fit]` to file `ifile.pdf` in the `iofolder` together with the +`systime` of creation. """ function plot_jold(jvals::Dict{Symbol,Any},systime::DateTime,iofolder::String,ifile::String) @@ -49,13 +50,11 @@ end #function plot_j """ - wrt_params(rxn, fit, sigma, rmse, R2, iofolder, time) + wrt_params(jvals, iofolder, systime) -For each reaction (rxn), from param in fit and sigma, print parameters -and 95% confidence together with RMSE and R^2 to file 'parameters.dat' -in the designated output folder (iofolder). - -Print the time of creation (time) to the output file. +Write the parameters, statistical data, and reaction labels stored in the +dictionary `jvals` to the file `parameters.dat` in the `iofolder` and state +the `systime` of creation. """ function wrt_params(jvals, iofolder, systime) diff --git a/src/rdfiles.jl b/src/rdfiles.jl index e313548..fa8d7a9 100644 --- a/src/rdfiles.jl +++ b/src/rdfiles.jl @@ -1,8 +1,11 @@ """ - get_fnames(scen) + get_fnames(scen, output) -From the scenario name (scen) used in TUV derive names for the TUV output file +From the scenario name (scen) used in TUV, derive names for the TUV output file and the output folder. + +If `output` is `true`, create a output folder `./params_`. Ask to overwrite, +if folder already exists. """ function setup_files(scen, output) # Test script argument and ask for input file, if missing