diff --git a/MiniObserve/.gitignore b/MiniObserve/.gitignore deleted file mode 100644 index 29126e4..0000000 --- a/MiniObserve/.gitignore +++ /dev/null @@ -1,24 +0,0 @@ -# Files generated by invoking Julia with --code-coverage -*.jl.cov -*.jl.*.cov - -# Files generated by invoking Julia with --track-allocation -*.jl.mem - -# System-specific files and directories generated by the BinaryProvider and BinDeps packages -# They contain absolute paths specific to the host computer, and so should not be committed -deps/deps.jl -deps/build.log -deps/downloads/ -deps/usr/ -deps/src/ - -# Build artifacts for creating documentation generated by the Documenter package -docs/build/ -docs/site/ - -# File generated by Pkg, the package manager, based on a corresponding Project.toml -# It records a fixed state of all packages used by the project. As such, it should not be -# committed for packages, but should be committed for applications that require a static -# environment. -Manifest.toml diff --git a/MiniObserve/LICENSE b/MiniObserve/LICENSE deleted file mode 100644 index f288702..0000000 --- a/MiniObserve/LICENSE +++ /dev/null @@ -1,674 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 3, 29 June 2007 - - Copyright (C) 2007 Free Software Foundation, Inc. - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The GNU General Public License is a free, copyleft license for -software and other kinds of works. - - The licenses for most software and other practical works are designed -to take away your freedom to share and change the works. By contrast, -the GNU General Public License is intended to guarantee your freedom to -share and change all versions of a program--to make sure it remains free -software for all its users. We, the Free Software Foundation, use the -GNU General Public License for most of our software; it applies also to -any other work released this way by its authors. You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -them if you wish), that you receive source code or can get it if you -want it, that you can change the software or use pieces of it in new -free programs, and that you know you can do these things. - - To protect your rights, we need to prevent others from denying you -these rights or asking you to surrender the rights. Therefore, you have -certain responsibilities if you distribute copies of the software, or if -you modify it: responsibilities to respect the freedom of others. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must pass on to the recipients the same -freedoms that you received. You must make sure that they, too, receive -or can get the source code. And you must show them these terms so they -know their rights. - - Developers that use the GNU GPL protect your rights with two steps: -(1) assert copyright on the software, and (2) offer you this License -giving you legal permission to copy, distribute and/or modify it. - - For the developers' and authors' protection, the GPL clearly explains -that there is no warranty for this free software. For both users' and -authors' sake, the GPL requires that modified versions be marked as -changed, so that their problems will not be attributed erroneously to -authors of previous versions. - - Some devices are designed to deny users access to install or run -modified versions of the software inside them, although the manufacturer -can do so. This is fundamentally incompatible with the aim of -protecting users' freedom to change the software. The systematic -pattern of such abuse occurs in the area of products for individuals to -use, which is precisely where it is most unacceptable. Therefore, we -have designed this version of the GPL to prohibit the practice for those -products. If such problems arise substantially in other domains, we -stand ready to extend this provision to those domains in future versions -of the GPL, as needed to protect the freedom of users. - - Finally, every program is threatened constantly by software patents. -States should not allow patents to restrict development and use of -software on general-purpose computers, but in those that do, we wish to -avoid the special danger that patents applied to a free program could -make it effectively proprietary. To prevent this, the GPL assures that -patents cannot be used to render the program non-free. - - The precise terms and conditions for copying, distribution and -modification follow. - - TERMS AND CONDITIONS - - 0. Definitions. - - "This License" refers to version 3 of the GNU General Public License. - - "Copyright" also means copyright-like laws that apply to other kinds of -works, such as semiconductor masks. - - "The Program" refers to any copyrightable work licensed under this -License. Each licensee is addressed as "you". "Licensees" and -"recipients" may be individuals or organizations. - - To "modify" a work means to copy from or adapt all or part of the work -in a fashion requiring copyright permission, other than the making of an -exact copy. The resulting work is called a "modified version" of the -earlier work or a work "based on" the earlier work. - - A "covered work" means either the unmodified Program or a work based -on the Program. - - To "propagate" a work means to do anything with it that, without -permission, would make you directly or secondarily liable for -infringement under applicable copyright law, except executing it on a -computer or modifying a private copy. Propagation includes copying, -distribution (with or without modification), making available to the -public, and in some countries other activities as well. - - To "convey" a work means any kind of propagation that enables other -parties to make or receive copies. Mere interaction with a user through -a computer network, with no transfer of a copy, is not conveying. - - An interactive user interface displays "Appropriate Legal Notices" -to the extent that it includes a convenient and prominently visible -feature that (1) displays an appropriate copyright notice, and (2) -tells the user that there is no warranty for the work (except to the -extent that warranties are provided), that licensees may convey the -work under this License, and how to view a copy of this License. If -the interface presents a list of user commands or options, such as a -menu, a prominent item in the list meets this criterion. - - 1. Source Code. - - The "source code" for a work means the preferred form of the work -for making modifications to it. "Object code" means any non-source -form of a work. - - A "Standard Interface" means an interface that either is an official -standard defined by a recognized standards body, or, in the case of -interfaces specified for a particular programming language, one that -is widely used among developers working in that language. - - The "System Libraries" of an executable work include anything, other -than the work as a whole, that (a) is included in the normal form of -packaging a Major Component, but which is not part of that Major -Component, and (b) serves only to enable use of the work with that -Major Component, or to implement a Standard Interface for which an -implementation is available to the public in source code form. A -"Major Component", in this context, means a major essential component -(kernel, window system, and so on) of the specific operating system -(if any) on which the executable work runs, or a compiler used to -produce the work, or an object code interpreter used to run it. - - The "Corresponding Source" for a work in object code form means all -the source code needed to generate, install, and (for an executable -work) run the object code and to modify the work, including scripts to -control those activities. However, it does not include the work's -System Libraries, or general-purpose tools or generally available free -programs which are used unmodified in performing those activities but -which are not part of the work. For example, Corresponding Source -includes interface definition files associated with source files for -the work, and the source code for shared libraries and dynamically -linked subprograms that the work is specifically designed to require, -such as by intimate data communication or control flow between those -subprograms and other parts of the work. - - The Corresponding Source need not include anything that users -can regenerate automatically from other parts of the Corresponding -Source. - - The Corresponding Source for a work in source code form is that -same work. - - 2. Basic Permissions. - - All rights granted under this License are granted for the term of -copyright on the Program, and are irrevocable provided the stated -conditions are met. This License explicitly affirms your unlimited -permission to run the unmodified Program. The output from running a -covered work is covered by this License only if the output, given its -content, constitutes a covered work. This License acknowledges your -rights of fair use or other equivalent, as provided by copyright law. - - You may make, run and propagate covered works that you do not -convey, without conditions so long as your license otherwise remains -in force. You may convey covered works to others for the sole purpose -of having them make modifications exclusively for you, or provide you -with facilities for running those works, provided that you comply with -the terms of this License in conveying all material for which you do -not control copyright. Those thus making or running the covered works -for you must do so exclusively on your behalf, under your direction -and control, on terms that prohibit them from making any copies of -your copyrighted material outside their relationship with you. - - Conveying under any other circumstances is permitted solely under -the conditions stated below. Sublicensing is not allowed; section 10 -makes it unnecessary. - - 3. Protecting Users' Legal Rights From Anti-Circumvention Law. - - No covered work shall be deemed part of an effective technological -measure under any applicable law fulfilling obligations under article -11 of the WIPO copyright treaty adopted on 20 December 1996, or -similar laws prohibiting or restricting circumvention of such -measures. - - When you convey a covered work, you waive any legal power to forbid -circumvention of technological measures to the extent such circumvention -is effected by exercising rights under this License with respect to -the covered work, and you disclaim any intention to limit operation or -modification of the work as a means of enforcing, against the work's -users, your or third parties' legal rights to forbid circumvention of -technological measures. - - 4. Conveying Verbatim Copies. - - You may convey verbatim copies of the Program's source code as you -receive it, in any medium, provided that you conspicuously and -appropriately publish on each copy an appropriate copyright notice; -keep intact all notices stating that this License and any -non-permissive terms added in accord with section 7 apply to the code; -keep intact all notices of the absence of any warranty; and give all -recipients a copy of this License along with the Program. - - You may charge any price or no price for each copy that you convey, -and you may offer support or warranty protection for a fee. - - 5. Conveying Modified Source Versions. - - You may convey a work based on the Program, or the modifications to -produce it from the Program, in the form of source code under the -terms of section 4, provided that you also meet all of these conditions: - - a) The work must carry prominent notices stating that you modified - it, and giving a relevant date. - - b) The work must carry prominent notices stating that it is - released under this License and any conditions added under section - 7. This requirement modifies the requirement in section 4 to - "keep intact all notices". - - c) You must license the entire work, as a whole, under this - License to anyone who comes into possession of a copy. This - License will therefore apply, along with any applicable section 7 - additional terms, to the whole of the work, and all its parts, - regardless of how they are packaged. This License gives no - permission to license the work in any other way, but it does not - invalidate such permission if you have separately received it. - - d) If the work has interactive user interfaces, each must display - Appropriate Legal Notices; however, if the Program has interactive - interfaces that do not display Appropriate Legal Notices, your - work need not make them do so. - - A compilation of a covered work with other separate and independent -works, which are not by their nature extensions of the covered work, -and which are not combined with it such as to form a larger program, -in or on a volume of a storage or distribution medium, is called an -"aggregate" if the compilation and its resulting copyright are not -used to limit the access or legal rights of the compilation's users -beyond what the individual works permit. Inclusion of a covered work -in an aggregate does not cause this License to apply to the other -parts of the aggregate. - - 6. Conveying Non-Source Forms. - - You may convey a covered work in object code form under the terms -of sections 4 and 5, provided that you also convey the -machine-readable Corresponding Source under the terms of this License, -in one of these ways: - - a) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by the - Corresponding Source fixed on a durable physical medium - customarily used for software interchange. - - b) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by a - written offer, valid for at least three years and valid for as - long as you offer spare parts or customer support for that product - model, to give anyone who possesses the object code either (1) a - copy of the Corresponding Source for all the software in the - product that is covered by this License, on a durable physical - medium customarily used for software interchange, for a price no - more than your reasonable cost of physically performing this - conveying of source, or (2) access to copy the - Corresponding Source from a network server at no charge. - - c) Convey individual copies of the object code with a copy of the - written offer to provide the Corresponding Source. This - alternative is allowed only occasionally and noncommercially, and - only if you received the object code with such an offer, in accord - with subsection 6b. - - d) Convey the object code by offering access from a designated - place (gratis or for a charge), and offer equivalent access to the - Corresponding Source in the same way through the same place at no - further charge. You need not require recipients to copy the - Corresponding Source along with the object code. If the place to - copy the object code is a network server, the Corresponding Source - may be on a different server (operated by you or a third party) - that supports equivalent copying facilities, provided you maintain - clear directions next to the object code saying where to find the - Corresponding Source. Regardless of what server hosts the - Corresponding Source, you remain obligated to ensure that it is - available for as long as needed to satisfy these requirements. - - e) Convey the object code using peer-to-peer transmission, provided - you inform other peers where the object code and Corresponding - Source of the work are being offered to the general public at no - charge under subsection 6d. - - A separable portion of the object code, whose source code is excluded -from the Corresponding Source as a System Library, need not be -included in conveying the object code work. - - A "User Product" is either (1) a "consumer product", which means any -tangible personal property which is normally used for personal, family, -or household purposes, or (2) anything designed or sold for incorporation -into a dwelling. In determining whether a product is a consumer product, -doubtful cases shall be resolved in favor of coverage. For a particular -product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status -of the particular user or of the way in which the particular user -actually uses, or expects or is expected to use, the product. A product -is a consumer product regardless of whether the product has substantial -commercial, industrial or non-consumer uses, unless such uses represent -the only significant mode of use of the product. - - "Installation Information" for a User Product means any methods, -procedures, authorization keys, or other information required to install -and execute modified versions of a covered work in that User Product from -a modified version of its Corresponding Source. The information must -suffice to ensure that the continued functioning of the modified object -code is in no case prevented or interfered with solely because -modification has been made. - - If you convey an object code work under this section in, or with, or -specifically for use in, a User Product, and the conveying occurs as -part of a transaction in which the right of possession and use of the -User Product is transferred to the recipient in perpetuity or for a -fixed term (regardless of how the transaction is characterized), the -Corresponding Source conveyed under this section must be accompanied -by the Installation Information. But this requirement does not apply -if neither you nor any third party retains the ability to install -modified object code on the User Product (for example, the work has -been installed in ROM). - - The requirement to provide Installation Information does not include a -requirement to continue to provide support service, warranty, or updates -for a work that has been modified or installed by the recipient, or for -the User Product in which it has been modified or installed. Access to a -network may be denied when the modification itself materially and -adversely affects the operation of the network or violates the rules and -protocols for communication across the network. - - Corresponding Source conveyed, and Installation Information provided, -in accord with this section must be in a format that is publicly -documented (and with an implementation available to the public in -source code form), and must require no special password or key for -unpacking, reading or copying. - - 7. Additional Terms. - - "Additional permissions" are terms that supplement the terms of this -License by making exceptions from one or more of its conditions. -Additional permissions that are applicable to the entire Program shall -be treated as though they were included in this License, to the extent -that they are valid under applicable law. If additional permissions -apply only to part of the Program, that part may be used separately -under those permissions, but the entire Program remains governed by -this License without regard to the additional permissions. - - When you convey a copy of a covered work, you may at your option -remove any additional permissions from that copy, or from any part of -it. (Additional permissions may be written to require their own -removal in certain cases when you modify the work.) You may place -additional permissions on material, added by you to a covered work, -for which you have or can give appropriate copyright permission. - - Notwithstanding any other provision of this License, for material you -add to a covered work, you may (if authorized by the copyright holders of -that material) supplement the terms of this License with terms: - - a) Disclaiming warranty or limiting liability differently from the - terms of sections 15 and 16 of this License; or - - b) Requiring preservation of specified reasonable legal notices or - author attributions in that material or in the Appropriate Legal - Notices displayed by works containing it; or - - c) Prohibiting misrepresentation of the origin of that material, or - requiring that modified versions of such material be marked in - reasonable ways as different from the original version; or - - d) Limiting the use for publicity purposes of names of licensors or - authors of the material; or - - e) Declining to grant rights under trademark law for use of some - trade names, trademarks, or service marks; or - - f) Requiring indemnification of licensors and authors of that - material by anyone who conveys the material (or modified versions of - it) with contractual assumptions of liability to the recipient, for - any liability that these contractual assumptions directly impose on - those licensors and authors. - - All other non-permissive additional terms are considered "further -restrictions" within the meaning of section 10. If the Program as you -received it, or any part of it, contains a notice stating that it is -governed by this License along with a term that is a further -restriction, you may remove that term. If a license document contains -a further restriction but permits relicensing or conveying under this -License, you may add to a covered work material governed by the terms -of that license document, provided that the further restriction does -not survive such relicensing or conveying. - - If you add terms to a covered work in accord with this section, you -must place, in the relevant source files, a statement of the -additional terms that apply to those files, or a notice indicating -where to find the applicable terms. - - Additional terms, permissive or non-permissive, may be stated in the -form of a separately written license, or stated as exceptions; -the above requirements apply either way. - - 8. Termination. - - You may not propagate or modify a covered work except as expressly -provided under this License. Any attempt otherwise to propagate or -modify it is void, and will automatically terminate your rights under -this License (including any patent licenses granted under the third -paragraph of section 11). - - However, if you cease all violation of this License, then your -license from a particular copyright holder is reinstated (a) -provisionally, unless and until the copyright holder explicitly and -finally terminates your license, and (b) permanently, if the copyright -holder fails to notify you of the violation by some reasonable means -prior to 60 days after the cessation. - - Moreover, your license from a particular copyright holder is -reinstated permanently if the copyright holder notifies you of the -violation by some reasonable means, this is the first time you have -received notice of violation of this License (for any work) from that -copyright holder, and you cure the violation prior to 30 days after -your receipt of the notice. - - Termination of your rights under this section does not terminate the -licenses of parties who have received copies or rights from you under -this License. If your rights have been terminated and not permanently -reinstated, you do not qualify to receive new licenses for the same -material under section 10. - - 9. Acceptance Not Required for Having Copies. - - You are not required to accept this License in order to receive or -run a copy of the Program. Ancillary propagation of a covered work -occurring solely as a consequence of using peer-to-peer transmission -to receive a copy likewise does not require acceptance. However, -nothing other than this License grants you permission to propagate or -modify any covered work. These actions infringe copyright if you do -not accept this License. Therefore, by modifying or propagating a -covered work, you indicate your acceptance of this License to do so. - - 10. Automatic Licensing of Downstream Recipients. - - Each time you convey a covered work, the recipient automatically -receives a license from the original licensors, to run, modify and -propagate that work, subject to this License. You are not responsible -for enforcing compliance by third parties with this License. - - An "entity transaction" is a transaction transferring control of an -organization, or substantially all assets of one, or subdividing an -organization, or merging organizations. If propagation of a covered -work results from an entity transaction, each party to that -transaction who receives a copy of the work also receives whatever -licenses to the work the party's predecessor in interest had or could -give under the previous paragraph, plus a right to possession of the -Corresponding Source of the work from the predecessor in interest, if -the predecessor has it or can get it with reasonable efforts. - - You may not impose any further restrictions on the exercise of the -rights granted or affirmed under this License. For example, you may -not impose a license fee, royalty, or other charge for exercise of -rights granted under this License, and you may not initiate litigation -(including a cross-claim or counterclaim in a lawsuit) alleging that -any patent claim is infringed by making, using, selling, offering for -sale, or importing the Program or any portion of it. - - 11. Patents. - - A "contributor" is a copyright holder who authorizes use under this -License of the Program or a work on which the Program is based. The -work thus licensed is called the contributor's "contributor version". - - A contributor's "essential patent claims" are all patent claims -owned or controlled by the contributor, whether already acquired or -hereafter acquired, that would be infringed by some manner, permitted -by this License, of making, using, or selling its contributor version, -but do not include claims that would be infringed only as a -consequence of further modification of the contributor version. For -purposes of this definition, "control" includes the right to grant -patent sublicenses in a manner consistent with the requirements of -this License. - - Each contributor grants you a non-exclusive, worldwide, royalty-free -patent license under the contributor's essential patent claims, to -make, use, sell, offer for sale, import and otherwise run, modify and -propagate the contents of its contributor version. - - In the following three paragraphs, a "patent license" is any express -agreement or commitment, however denominated, not to enforce a patent -(such as an express permission to practice a patent or covenant not to -sue for patent infringement). To "grant" such a patent license to a -party means to make such an agreement or commitment not to enforce a -patent against the party. - - If you convey a covered work, knowingly relying on a patent license, -and the Corresponding Source of the work is not available for anyone -to copy, free of charge and under the terms of this License, through a -publicly available network server or other readily accessible means, -then you must either (1) cause the Corresponding Source to be so -available, or (2) arrange to deprive yourself of the benefit of the -patent license for this particular work, or (3) arrange, in a manner -consistent with the requirements of this License, to extend the patent -license to downstream recipients. "Knowingly relying" means you have -actual knowledge that, but for the patent license, your conveying the -covered work in a country, or your recipient's use of the covered work -in a country, would infringe one or more identifiable patents in that -country that you have reason to believe are valid. - - If, pursuant to or in connection with a single transaction or -arrangement, you convey, or propagate by procuring conveyance of, a -covered work, and grant a patent license to some of the parties -receiving the covered work authorizing them to use, propagate, modify -or convey a specific copy of the covered work, then the patent license -you grant is automatically extended to all recipients of the covered -work and works based on it. - - A patent license is "discriminatory" if it does not include within -the scope of its coverage, prohibits the exercise of, or is -conditioned on the non-exercise of one or more of the rights that are -specifically granted under this License. You may not convey a covered -work if you are a party to an arrangement with a third party that is -in the business of distributing software, under which you make payment -to the third party based on the extent of your activity of conveying -the work, and under which the third party grants, to any of the -parties who would receive the covered work from you, a discriminatory -patent license (a) in connection with copies of the covered work -conveyed by you (or copies made from those copies), or (b) primarily -for and in connection with specific products or compilations that -contain the covered work, unless you entered into that arrangement, -or that patent license was granted, prior to 28 March 2007. - - Nothing in this License shall be construed as excluding or limiting -any implied license or other defenses to infringement that may -otherwise be available to you under applicable patent law. - - 12. No Surrender of Others' Freedom. - - If conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot convey a -covered work so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you may -not convey it at all. For example, if you agree to terms that obligate you -to collect a royalty for further conveying from those to whom you convey -the Program, the only way you could satisfy both those terms and this -License would be to refrain entirely from conveying the Program. - - 13. Use with the GNU Affero General Public License. - - Notwithstanding any other provision of this License, you have -permission to link or combine any covered work with a work licensed -under version 3 of the GNU Affero General Public License into a single -combined work, and to convey the resulting work. The terms of this -License will continue to apply to the part which is the covered work, -but the special requirements of the GNU Affero General Public License, -section 13, concerning interaction through a network will apply to the -combination as such. - - 14. Revised Versions of this License. - - The Free Software Foundation may publish revised and/or new versions of -the GNU General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - - Each version is given a distinguishing version number. If the -Program specifies that a certain numbered version of the GNU General -Public License "or any later version" applies to it, you have the -option of following the terms and conditions either of that numbered -version or of any later version published by the Free Software -Foundation. If the Program does not specify a version number of the -GNU General Public License, you may choose any version ever published -by the Free Software Foundation. - - If the Program specifies that a proxy can decide which future -versions of the GNU General Public License can be used, that proxy's -public statement of acceptance of a version permanently authorizes you -to choose that version for the Program. - - Later license versions may give you additional or different -permissions. However, no additional obligations are imposed on any -author or copyright holder as a result of your choosing to follow a -later version. - - 15. Disclaimer of Warranty. - - THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY -APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT -HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY -OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, -THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM -IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF -ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. Limitation of Liability. - - IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS -THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY -GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE -USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF -DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD -PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), -EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF -SUCH DAMAGES. - - 17. Interpretation of Sections 15 and 16. - - If the disclaimer of warranty and limitation of liability provided -above cannot be given local legal effect according to their terms, -reviewing courts shall apply local law that most closely approximates -an absolute waiver of all civil liability in connection with the -Program, unless a warranty or assumption of liability accompanies a -copy of the Program in return for a fee. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -state the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - - Copyright (C) - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . - -Also add information on how to contact you by electronic and paper mail. - - If the program does terminal interaction, make it output a short -notice like this when it starts in an interactive mode: - - Copyright (C) - This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. - -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, your program's commands -might be different; for a GUI interface, you would use an "about box". - - You should also get your employer (if you work as a programmer) or school, -if any, to sign a "copyright disclaimer" for the program, if necessary. -For more information on this, and how to apply and follow the GNU GPL, see -. - - The GNU General Public License does not permit incorporating your program -into proprietary programs. If your program is a subroutine library, you -may consider it more useful to permit linking proprietary applications with -the library. If this is what you want to do, use the GNU Lesser General -Public License instead of this License. But first, please read -. diff --git a/MiniObserve/README.md b/MiniObserve/README.md deleted file mode 100644 index b734571..0000000 --- a/MiniObserve/README.md +++ /dev/null @@ -1,40 +0,0 @@ -# MiniObserve.jl -Minimalist (and minimally intrusive) macro set for extracting information from complex objects. - - - -@observe(statstype, model [, user_arg1...], declarations) - -Generate a full analysis suite for a model. - -Given a declaration - -```Julia -@observe Data model user1 user2 begin - @record "time" model.time - @record "N" Int length(model.population) - - @for ind in model.population begin - @stat("capital", MaxMinAcc{Float64}, MeanVarAcc{FloatT}) <| ind.capital - @stat("n_alone", CountAcc) <| has_neighbours(ind) - end - - @record u1 user1 - @record u2 user1 * user2 -end -``` - -a type Data will be generated that provides (at least) the following members: - -```Julia -struct Data - time :: Float64 - N :: Int - capital :: @NamedTuple{max :: Float64, min :: Float64, mean :: Float64, var :: Float64} - n_alone :: @NamedTuple{N :: Int} - u1 :: Float64 - u2 :: Float64 -end -``` - -The macro will also create a method for `observe(::Type{Data), model...)` that will perform the required calculations and returns a `Data` object. Use `print_header` to print a header for the generated type to an output and `log_results` to print the content of a data object. diff --git a/MiniObserve/src/MiniObserve.jl b/MiniObserve/src/MiniObserve.jl deleted file mode 100644 index 73c9a39..0000000 --- a/MiniObserve/src/MiniObserve.jl +++ /dev/null @@ -1,28 +0,0 @@ -module MiniObserve - -using Reexport - -include("Observation.jl") -include("StatsAccumulator.jl") - -@reexport using .Observation -@reexport using .StatsAccumulator - - -# these are not needed any more - -#import .Observation.prefixes -#prefixes(::Type{CountAcc}) = ["#"] -#prefixes(::Type{<:MeanVarAcc}) = ["mean", "var"] -#prefixes(::Type{<:MaxMinAcc}) = ["max", "min"] -# - -#import .Observation.printable -#printable(acc :: CountAcc, FS="\t") = acc.n -#printable(acc :: MaxMinAcc, FS="\t") = (acc.max, FS, acc.min) -#function printable(acc :: MeanVarAcc, FS="\t") -# res = result(acc) -# (res[1], FS, res[2]) -#end - -end # module diff --git a/MiniObserve/src/Observation.jl b/MiniObserve/src/Observation.jl deleted file mode 100644 index 2c494e2..0000000 --- a/MiniObserve/src/Observation.jl +++ /dev/null @@ -1,330 +0,0 @@ -module Observation - -export print_header, @observe, log_results - -using MacroTools - -"obtain a named tuple type with the same field types and names as `struct_T`" -tuple_type(struct_T) = NamedTuple{fieldnames(struct_T), Tuple{fieldtypes(struct_T)...}} - -"construct a named tuple from `x`" -@generated function to_named_tuple(x) - if x <: NamedTuple - return :x - end - - # constructor call - tuptyp = Expr(:quote, tuple_type(x)) - - # constructor arguments - tup = Expr(:tuple) - for i in 1:fieldcount(x) - push!(tup.args, :(getfield(x, $i)) ) - end - - # both put together - :($tuptyp($tup)) -end - - -"translate accumulator types into prefixes for the header (e.g. min, max, etc.)" -stat_names(::Type{T}) where {T} = fieldnames(result_type(T)) - -# We could make this a generated function as well, but since the header -# is only printed once at the beginning, the additional time needed for -# runtime introspection is worth the reduction in complexity. -"Print a header for an observation type `stats_t` to `output` using field separator `FS`, name separator `NS` and line separator `LS`." -function print_header(output, stats_t; FS="\t", NS="_", LS="\n") - fn = fieldnames(stats_t) - ft = fieldtypes(stats_t) - - for (i, (name, typ)) in enumerate(zip(fn, ft)) - if typ <: NamedTuple - # aggregate stat - header(output, string(name), string.(fieldnames(typ)), FS, NS) - else - # single stat - print(output, string(name)) - end - - if i < length(fn) - print(output, FS) - end - end - - print(output, LS) -end - -# print header for aggregate stat -function header(out, stat_name, stat_names, FS, NS) - @assert length(stat_names) > 0 - - print(out, join((stat_name * NS) .* stat_names, FS)) -end - -# It's quite possibly overkill to make this a generated function, but we -# don't want anybody accusing us of wasting CPU cycles. -"Print results stored in `stats` to `output` using field separator `FS` and line separator `LS`." -@generated function log_results(out, stats; FS="\t", LS="\n") - fn = fieldnames(stats) - ft = fieldtypes(stats) - - fn_body = Expr(:block) - - # all fields of stats - for (i, (name, typ)) in enumerate(zip(fn, ft)) - # aggregate stats - if typ <: NamedTuple - # go through all elements of stats.name - for (j, tname) in enumerate(fieldnames(typ)) - push!(fn_body.args, :(print(out, stats.$name.$tname))) - if j < length(fieldnames(typ)) - push!(fn_body.args, :(print(out, FS))) - end - end - # single values - else - push!(fn_body.args, :(print(out, stats.$name))) - end - - if i < length(fn) - push!(fn_body.args, :(print(out, FS))) - end - end - - push!(fn_body.args, :(print(out, LS))) - - fn_body -end - - -# process a single stats declaration (@record) -function process_single(name, typ, expr) - # no type specified, default to float - if typ == nothing - typ = :Float64 - end - - tmp_name = gensym("tmp") - :($tmp_name = $(esc(expr))) - - [:($name :: $(esc(typ)))], # type - [:($tmp_name = $(esc(expr)))], # body - [:($tmp_name)] # constructor -end - -# it would be much nicer to use a generated function for this, but -# unfortunately we are already operating on types -"concatenate named tuple and/or struct types into one single named tuple" -function joined_named_tuple_T(types...) - ns = Expr(:tuple) - ts = Expr(:curly) - push!(ts.args, :Tuple) - - for t in types - fnames = fieldnames(t) - ftypes = fieldtypes(t) - - append!(ns.args, QuoteNode.(fnames)) - append!(ts.args, ftypes) - end - - ret = :(NamedTuple{}) - push!(ret.args, ns) - push!(ret.args, ts) - - eval(ret) -end - - -# process an aggregate stats declaration (@for) -function process_aggregate(var, collection, decls) - stat_type_code = [] - body_code = [] - res_code = [] - - decl_code = [] - loop_code = [] - - lines = rmlines(decls).args - - for (i, line) in enumerate(lines) - @capture(line, @stat(statname_String, stattypes__) <| expr_) || - error("expected:@stat(, {, }) <| ") - - # code to declare stat property in stats struct - # creates a single named tuple type from all result types of all stats - sname = Symbol(statname) - prop_code = :($sname :: joined_named_tuple_T()) - for t in stattypes - # not a constructor call - if ! @capture(t, typ_(args__)) - typ = t - end - - push!(prop_code.args[2].args, :($(esc(:result_type))($(esc(typ))))) - end - - push!(stat_type_code, prop_code) - - # code to store result of user code (to be fed into stats objects) - # (inside loop) - tmp_name = gensym("tmp") - push!(loop_code, :($tmp_name = $(esc(expr)))) - - # expression that merges all results for this stat into single named tuple - res_expr = length(stattypes) > 1 ? :(merge()) : :(identity()) - - # all stats for this specific term - for (j, stattype) in enumerate(stattypes) - # declaration of accumulator - vname = gensym("stat") - # goes into main body (outside of loop) - if @capture(stattype, typ_(args__)) - # paste in constructor expression - push!(body_code, :($(esc(vname)) = $(esc(stattype)))) - else - # create constructor call - push!(body_code, :($(esc(vname)) = $(esc(stattype))())) - end - - # add value to accumulator - push!(loop_code, :($(esc(:add!))($(esc(vname)), $tmp_name))) - # add to named tuple argument of constructor call - push!(res_expr.args, :(to_named_tuple($(esc(:results))($(esc(vname)))))) - end - - # another argument for the main constructor call - push!(res_code, res_expr) - - end - - # add the loop to the main body - push!(body_code, :(for $(esc(var)) in $(esc(collection)); $(loop_code...); end)) - - stat_type_code, body_code, res_code -end - - -""" - -@observe(statstype, model [, user_arg1...], declarations) - -Generate a full analysis suite for a model. - -Given a declaration - -```Julia -@observe Data model user1 user2 begin - @record "time" model.time - @record "N" Int length(model.population) - - @for ind in model.population begin - @stat("capital", MaxMinAcc{Float64}, MeanVarAcc{FloatT}) <| ind.capital - @stat("n_alone", CountAcc) <| has_neighbours(ind) - end - - @record u1 user1 - @record u2 user1 * user2 -end -``` - -a type Data will be generated that provides (at least) the following members: - -```Julia -struct Data - time :: Float64 - N :: Int - capital :: @NamedTuple{max :: Float64, min :: Float64, mean :: Float64, var :: Float64} - n_alone :: @NamedTuple{N :: Int} - u1 :: Float64 - u2 :: Float64 -end -``` - -The macro will also create a method for `observe(::Type{Data), model...)` that will perform the required calculations and returns a `Data` object. Use `print_header` to print a header for the generated type to an output and `log_results` to print the content of a data object. -""" -macro observe(tname, model, args_and_decl...) - observe_syntax = "@observe [ ...] " - - if typeof(tname) != Symbol - error("usage: $observe_syntax") - end - - if length(args_and_decl) < 1 - error("usage: $observe_syntax") - end - - decl = args_and_decl[end] - - if typeof(decl) != Expr || decl.head != :block - error("usage: $observe_syntax") - end - - ana_func = if length(args_and_decl) == 1 - :(function $(esc(:observe))(::$(esc(:Type)){$(esc(tname))}, $(esc(model))); end) - else - :(function $(esc(:observe))(::$(esc(:Type)){$(esc(tname))}, $(esc(model)), - $(esc(args_and_decl[1:end-1]...))); end) - end - - ana_body = ana_func.args[2].args - - stats_type = :(struct $(esc(tname)); end) - - stats_constr = :($(esc(tname))()) - - syntax = "single or population stats declaration expected:\n" * - "\t@for in |" * - "\t@record |" * - "\t@record " - - # go through declaration expression by expression - # each expression is translated into three bits of code: - # * additional fields for the stats type - # * additional code to run during the analysis - # * additional arguments for the stats object constructor call - lines = rmlines(decl).args - for (i, line) in enumerate(lines) - if typeof(line) != Expr || line.head != :macrocall - dump(line) - error(syntax) - end - - # single stat - if line.args[1] == Symbol("@record") - typ = nothing - @capture(line, @record(name_String, expr_)) || - @capture(line, @record(name_String, typ_, expr_)) || - error("expecting: @record [] ") - stats_type_c, ana_body_c, stats_constr_c = - process_single(Symbol(name), typ, expr) - - # aggregate stat - elseif line.args[1] == Symbol("@for") - @capture(line, @for var_Symbol in expr_ begin block_ end) || - error("expecting: @for in ") - stats_type_c, ana_body_c, stats_constr_c = - process_aggregate(var, expr, block) - - else - error(syntax) - end - - # add code to respective bits - append!(ana_func.args[2].args, ana_body_c) - append!(stats_type.args[3].args, stats_type_c) - append!(stats_constr.args, stats_constr_c) - end - - # add constructor call as last line of analysis function - push!(ana_func.args[2].args, stats_constr) - - ret = Expr(:block) - push!(ret.args, stats_type) - push!(ret.args, ana_func) - - ret -end - -end # module diff --git a/MiniObserve/src/StatsAccumulator.jl b/MiniObserve/src/StatsAccumulator.jl deleted file mode 100644 index 0a5e2d3..0000000 --- a/MiniObserve/src/StatsAccumulator.jl +++ /dev/null @@ -1,130 +0,0 @@ -module StatsAccumulator - -export CountAcc, MeanVarAcc, add!, results, result_type, MaxMinAcc, HistAcc - - - -# per default the struct itself is assumed to contain the results (see e.g. min/max) - -"type of results of accumulator `T`, overload as needed" -result_type(::Type{T}) where {T} = T - -"results, overload as needed" -results(t :: T) where {T} = t - - -### CountAcc - -mutable struct CountAcc - n :: Int -end - -CountAcc() = CountAcc(0) - -add!(acc :: CountAcc, cond) = cond ? acc.n += 1 : 0 - - -mutable struct MeanVarAcc{T} - sum :: T - sum_sq :: T - n :: Int -end - -MeanVarAcc{T}() where {T} = MeanVarAcc(T(0), T(0), 0) - -function add!(acc :: MeanVarAcc{T}, v :: T) where {T} - acc.sum += v - acc.sum_sq += v*v - acc.n += 1 -end - - -results(acc :: MeanVarAcc{T}) where {T} = - (mean = acc.sum / acc.n, var = (acc.sum_sq - acc.sum*acc.sum/acc.n) / (acc.n - 1)) - -result_type(::Type{MeanVarAcc{T}}) where {T} = @NamedTuple{mean::T, var::T} - - - -#mutable struct MeanVarAcc2{T} -# m :: T -# m2 :: T -# n :: Int -#end -# -#MeanVarAcc2{T}() where T = MeanVarAcc2(T(0), T(0), 0) -# -#function add!(acc :: MeanVarAcc2{T}, v :: T) where T -# delta = v - acc.m -# acc.n += 1 -# delta_n = delta / acc.n -# acc.m += delta_n -# acc.m2 += delta * (delta - delta_n) -#end -# -#result(acc :: MeanVarAcc2{T}) where {T} = acc.m, acc.m2 / acc.n - - -mutable struct MaxMinAcc{T} - max :: T - min :: T -end - - -MaxMinAcc{T}() where {T} = MaxMinAcc(typemin(T), typemax(T)) - - -function add!(acc :: MaxMinAcc{T}, v :: T) where {T} - acc.max = max(acc.max, v) - acc.min = min(acc.min, v) -end - - -mutable struct HistAcc{T} - bins :: Vector{Int} - min :: T - width :: T -end - -HistAcc(min::T = T(0), width::T = T(1)) where {T} = HistAcc{T}([], min, width) - -function add!(acc :: HistAcc{T}, v :: T) where {T} - if v < acc.min - return acc - end - - bin = floor(Int, (v - acc.min) / acc.width) + 1 - n = length(acc.bins) - if bin > n - sizehint!(acc.bins, bin) - for i in (n+1):bin - push!(acc.bins, 0) - end - end - - acc.bins[bin] += 1 - - acc -end - -#results(acc::HistAcc{T}) where {T} = (;bins = acc.bins) -results(acc::HistAcc) = (;bins = acc.bins) - -#result_type(::Type{HistAcc{T}}) where {T} = @NamedTuple{bins::Vector{Int}} -result_type(::Type{HistAcc}) = @NamedTuple{bins::Vector{Int}} - -# does not work with results/result_type, maybe rework as tuples? - -#struct AccList -# list :: Vector{Any} -#end -# -#AccList() = AccList([]) -# -#function add!(al :: AccList, v :: T) where {T} -# for a in al.list -# add!(a, v) -# end -#end - -end # module diff --git a/MiniObserve/test/runtests.jl b/MiniObserve/test/runtests.jl deleted file mode 100644 index 6c81706..0000000 --- a/MiniObserve/test/runtests.jl +++ /dev/null @@ -1,95 +0,0 @@ -include("../src/MiniObserve.jl") -using .MiniObserve - -using Test - -struct Agent - capital :: Float64 - n :: Int -end - -has_neighbours(a) = a.n > 0 - -struct Model - time :: Float64 - population :: Vector{Agent} -end - -function Model() - t = 10.0 - p = [Agent(0.0, 2), Agent(2.0, 0), Agent(1.0, 10)] - - Model(t, p) -end - - -@observe Data model begin - @record "time" model.time - @record "N" Int length(model.population) - - @for ind in model.population begin - @stat("capital", MaxMinAcc{Float64}) <| ind.capital - @stat("n_alone", CountAcc) <| !has_neighbours(ind) - end -end - - -@testset "Analysis" begin - -m = Model() -result = observe(Data, m) - -@test typeof(result) == Data -@test result.time == m.time -@test result.N == length(m.population) -@test result.capital.max == 2.0 -@test result.capital.min == 0.0 -@test result.n_alone.n == 1 - -end - - -@testset "Output" begin - -m = Model() -result = observe(Data, m) - -io = IOBuffer() - -print_header(io, Data) -@test String(take!(io)) == "time\tN\tcapital_max\tcapital_min\tn_alone_n\n" - -log_results(io, result) -@test String(take!(io)) == "10.0\t3\t2.0\t0.0\t1\n" - -end - - -@observe Data2 model user1 begin - @record "time" model.time - @record "N" Int length(model.population) - - @for ind in model.population begin - @stat("capital", MaxMinAcc{Float64}) <| ind.capital - @stat("n_alone", CountAcc) <| !has_neighbours(ind) - end - - @record "user" user1 -end - - -@testset "User data" begin - -m = Model() -result = observe(Data2, m, 42) - -@test typeof(result) == Data2 -@test result.time == m.time -@test result.N == length(m.population) -@test result.capital.max == 2.0 -@test result.capital.min == 0.0 -@test result.n_alone.n == 1 -@test result.user == 42.0 - -end - diff --git a/analysis.jl b/analysis.jl index 4bd4c08..1d4ed87 100644 --- a/analysis.jl +++ b/analysis.jl @@ -5,23 +5,98 @@ const MVA = MeanVarAcc{Float64} # maximum, minimum const MMA = MaxMinAcc{Float64} -@observe Data model begin +function income_deciles(pop) + incomes = [ income(p) for p in pop ] + sort!(incomes) + + dec_size = length(pop) ÷ 10 + inc_decs = zeros(10) + + for i in 1:(10*dec_size) + inc = incomes[i] + inc_decs[(i-1) ÷ dec_size + 1] += inc + end + + inc_decs ./ dec_size +end + +@observe Data model t pars begin + @record "time" Float64 t + @record "N" Int length(model.pop) + @for house in model.houses begin # format: # @stat(name, accumulators...) <| expression - @stat("hh_size", MVA, MMA) <| Float64(length(house.occupants)) + @stat("hh_size", MVA, MMA, HistAcc(1.0, 1.0)) <| Float64(length(house.occupants)) + end + + # households with children + @for house in Iterators.filter(h->!isEmpty(h), model.houses) begin + i_c = findfirst(p->age(p)<18, house.occupants) + if i_c != nothing + child = house.occupants[i_c] + end + + is_lp = i_c != nothing && !isOrphan(child) && isSingle(guardians(child)[1]) + + # all hh with children + @stat("n_ch_hh", CountAcc) <| (i_c != nothing) + # hh with children with lone parents + @stat("n_lp_hh", CountAcc) <| is_lp + # number of children in lp households + @if is_lp @stat("n_ch_lp_hh", HistAcc(0, 1)) <| count(p->age(p)<18, house.occupants) + + ncs = netCareSupply(house) + scn = householdSocialCareNeed(house, model, pars.carepars) + cbn = careBalance(house) + unmet_pre = min(0, max(ncs, -scn)) + unmet_post = min(0, max(cbn, -scn)) + + @stat("care_supply", MVA) <| Float64(ncs) + @stat("unmet_care", MVA) <| Float64(min(cbn, 0)) + @stat("unmet_scare_pre", MVA) <| Float64(unmet_pre) + @stat("unmet_scare_post", MVA) <| Float64(unmet_post) end @for person in model.pop begin - @stat("married", CountAcc) <| (alive(person) && !isSingle(person)) - @stat("age", MVA) <| Float64(age(person)) - @stat("alive", CountAcc) <| alive(person) - @stat("eligible", CountAcc) <| (alive(person) && isFemale(person) && age(person) > 17) - @stat("eligible2", CountAcc) <| (alive(person) && isSingle(person) && isFemale(person) && age(person) > 17) + @stat("age", MVA, HistAcc(0.0, 3.0)) <| Float64(age(person)) + @stat("married", CountAcc) <| (!isSingle(person)) + @stat("single", CountAcc) <| (age(person) > 18 && isSingle(person)) + @stat("alive", CountAcc) <| true + @stat("class", HistAcc(0.0, 1.0, 4.0)) <| Float64(classRank(person)) + @stat("careneed", HistAcc(0.0, 1.0)) <| Float64(careNeedLevel(person)) + @stat("income", MVA) <| income(person) + @stat("n_orphans", CountAcc) <| isOrphan(person) + @if isFemale(person) @stat("f_status", HistAcc(0, 1, 5)) <| Int(status(person)) + @if isMale(person) @stat("m_status", HistAcc(0, 1, 5)) <| Int(status(person)) + @stat("p_care_supply", HistAcc(0.0, 4.0), MVA) <| + Float64(socialCareSupply(person, pars.carepars)) + @stat("p_care_demand", HistAcc(0.0, 4.0), MVA) <| + Float64(socialCareDemand(person, pars.carepars)) + end + + @for person in Iterators.filter(p -> isFemale(p) && ! isSingle(p), model.pop) begin + @stat("age_diff", HistAcc(-10.0, 1.0)) <| Float64(age(partner(person)) - age(person)) end + @record "income_deciles" Vector{Float64} income_deciles(model.pop) +end + + +function saveHouses(io, houses) + dump_header(io, houses[1]) + println(io) + for h in houses + Utilities.dump(io, h) + println(io) + end +end - @for person in Iterators.filter(p->alive(p), model.pop) begin - @stat("hist_age", HistAcc(0.0, 1.0)) <| Float64(age(person)) +function saveAgents(io, pop) + dump_header(io, pop[1]) + println(io) + for p in pop + Utilities.dump(io, p) + println(io) end end diff --git a/data/age_pyramid_1921_5y.csv b/data/age_pyramid_1921_5y.csv new file mode 100644 index 0000000..32dc3a7 --- /dev/null +++ b/data/age_pyramid_1921_5y.csv @@ -0,0 +1,20 @@ +0.093024440783328 0.0827936961738307 +0.097733700782601 0.0884520878319922 +0.101637660226789 0.0920024066878463 +0.0955906032556471 0.089606268291181 +0.080130890662082 0.0859637300834971 +0.0741323531047086 0.0817854918314955 +0.0708881359743016 0.0767055532504924 +0.0704455968742654 0.0742960387573657 +0.0676646101332325 0.0695618091750936 +0.0642955813751619 0.0627903243880057 +0.0537210600645447 0.0526528584970517 +0.0432419178523725 0.0428598901847718 +0.0332629073397038 0.0343623337199782 +0.0248606947880468 0.0270903305460577 +0.0155179690846688 0.0189950664918184 +0.00877111500434379 0.0118132636362994 +0.00370656233093239 0.00569655138995309 +0.00112479840515525 0.00201938675897687 +0.000216152051986698 0.000467507190282796 +3.32499061284888E-05 8.54051140097701E-05 diff --git a/data/birthrate_early.csv b/data/birthrate_early.csv new file mode 100644 index 0000000..e9f432a --- /dev/null +++ b/data/birthrate_early.csv @@ -0,0 +1,74 @@ +1887 0.0284689240992265 +1888 0.0280397925235824 +1889 0.0279504234928952 +1890 0.0272473375875021 +1891 0.0283416396842529 +1892 0.0276184569715819 +1893 0.0278774392509282 +1894 0.0269082660174837 +1895 0.0275496352728506 +1896 0.0271928402413186 +1897 0.0270808231613537 +1898 0.0268992194310167 +1899 0.0267596500736984 +1900 0.026473030961302 +1901 0.0263078563827994 +1902 0.0263406989761939 +1903 0.026347351029432 +1904 0.0260386187733799 +1905 0.0254092059710382 +1906 0.0253331903480308 +1907 0.0246434662923147 +1908 0.0249830023706027 +1909 0.0241193409629488 +1910 0.0234046295409866 +1911 0.0244939535148306 +1912 0.0242091302131516 +1913 0.024242138165388 +1914 0.0240411668439762 +1915 0.0231344849679049 +1916 0.0227471426844581 +1917 0.0198773777233684 +1918 0.0198935627305341 +1919 0.0192390107139281 +1920 0.0258176848491514 +1921 0.0228162062313087 +1922 0.0208596596099839 +1923 0.0201989987231588 +1924 0.0192785436754772 +1925 0.0187034857904085 +1926 0.0182489561586639 +1927 0.0171139285085402 +1928 0.0171638369039634 +1929 0.0166618486365922 +1930 0.0167630734229846 +1931 0.0162777382275316 +1932 0.0157565339376281 +1933 0.014865864144454 +1934 0.0152539964856641 +1935 0.0151791928480749 +1936 0.0152954357674916 +1937 0.0153055704757595 +1938 0.0154876711001998 +1939 0.0152821692742236 +1940 0.0152484727086284 +1941 0.0155052328483811 +1942 0.0174142318886357 +1943 0.0167943872561208 +1944 0.0181243344971894 +1945 0.0163527016225968 +1946 0.0195000796116584 +1947 0.0206995137135209 +1948 0.018091627159566 +1949 0.0169934632731319 +1950 0.0162444746583567 +1951 0.0158419986119645 +1952 0.0157233705868822 +1953 0.0158968748579346 +1954 0.0156558764027901 +1955 0.0154931388271134 +1956 0.0161211523244796 +1957 0.0165557590676295 +1958 0.01685295000242 +1959 0.0169096144259695 +1960 0.0175337438541219 diff --git a/data/deathrates_early.csv b/data/deathrates_early.csv new file mode 100644 index 0000000..21e57d6 --- /dev/null +++ b/data/deathrates_early.csv @@ -0,0 +1,61 @@ +1900 0.0169086080292361 150 +1901 0.0157841697521799 147.50 +1902 0.0151971584548143 129.90 +1903 0.0145272282266502 129.10 +1904 0.0152846656059177 141.70 +1905 0.0143672523545397 126.10 +1906 0.0145281139085494 129.60 +1907 0.0142958950838863 116.10 +1908 0.0140837144579569 120.20 +1909 0.0138121497321398 108.40 +1910 0.0128705202389354 105.80 +1911 0.014716068812841 127.10 +1912 0.0137108246644137 96.20 +1913 0.0141033715886648 108.50 +1914 0.0142461397354092 105.30 +1915 0.0161097155567375 111.60 +1916 0.0147921985973066 91.80 +1917 0.0148166321859595 97.90 +1918 0.0180699813046334 97.60 +1919 0.0140226014749407 90.90 +1920 0.0127232944756039 81.80 +1921 0.0123938311000567 83.80 +1922 0.0130715179951051 80.10 +1923 0.0118227412365836 70.80 +1924 0.0125628486641596 78.20 +1925 0.012391918294849 77.30 +1926 0.0118628808251654 72.30 +1927 0.0125166182426484 72.40 +1928 0.0119166546111323 68 +1929 0.0136281953160876 76.30 +1930 0.0116991254965798 63.10 +1931 0.0124563307403806 68.50 +1932 0.0122582497032481 68.30 +1933 0.0124562983662941 66.40 +1934 0.0119588565550936 61.40 +1935 0.0119765727514215 60.40 +1936 0.0123391240258871 62.10 +1937 0.0126414823022885 61.10 +1938 0.0117824740336168 55.50 +1939 0.0122373321948275 53.60 +1940 0.0146266500395401 61 +1941 0.0139820014976466 63.30 +1942 0.0126876790830946 52.90 +1943 0.0121334974389577 51.90 +1944 0.0118360448703662 47.60 +1945 0.0116507050703838 48.80 +1946 0.0117041589946885 42.70 +1947 0.0121264385218021 43.70 +1948 0.0109127939048472 36 +1949 0.0117199340366772 34.10 +1950 0.0117133471611603 31.20 +1951 0.0125835157864175 31.10 +1952 0.0113784474074544 28.80 +1953 0.0114091107645539 27.60 +1954 0.0113936991897945 26.40 +1955 0.0116969895634798 25.80 +1956 0.0116830814617992 24.40 +1957 0.0114951915411568 23.90 +1958 0.0116943032767049 23.30 +1959 0.0116658615028399 23.10 +1960 0.0115199388992315 22.50 diff --git a/src/DeclUtils.jl b/lib/DeclUtils.jl similarity index 78% rename from src/DeclUtils.jl rename to lib/DeclUtils.jl index 20b1bdc..af9db65 100644 --- a/src/DeclUtils.jl +++ b/lib/DeclUtils.jl @@ -6,12 +6,12 @@ module DeclUtils function setter(field, type) name = setter_name(field) - :($(esc(name))(x::$type, value) = (x.$field = value)) + :($(esc(name))(x::$(esc(type)), value) = (x.$field = value)) end function getter(field, type) name = getter_name(field) - :($(esc(name))(x::$type) = x.$field) + :($(esc(name))(x::$(esc(type))) = x.$field) end macro decl_setters(type, fields...) @@ -46,6 +46,14 @@ module DeclUtils :(export $(esc(gname)), $(esc(sname))) end + macro export_forward(type, fieldvec) + @assert fieldvec.head == :vect + fields = fieldvec.args + Expr(:block, + [setter(f, type) for f in fields]..., + [getter(f, type) for f in fields]..., + [getset_export(f) for f in fields]...) + end macro export_forward(type, sub, fieldvec) @assert fieldvec.head == :vect diff --git a/src/ParamUtils.jl b/lib/ParamUtils.jl similarity index 100% rename from src/ParamUtils.jl rename to lib/ParamUtils.jl diff --git a/loadLibsPath.jl b/lib/loadLibsPath.jl similarity index 100% rename from loadLibsPath.jl rename to lib/loadLibsPath.jl diff --git a/lpm.jl b/lpm.jl deleted file mode 100644 index dad034c..0000000 --- a/lpm.jl +++ /dev/null @@ -1,5 +0,0 @@ -include("./loadLibsPath.jl") - -addToLoadPath!("src/generic") - -include("mainHelpers.jl") diff --git a/main.jl b/main.jl index 9adfc77..557d149 100644 --- a/main.jl +++ b/main.jl @@ -2,14 +2,13 @@ using Random # load main simulation code -include("lpm.jl") +include("mainHelpers.jl") # create parameters const simPars, pars = loadParameters(ARGS) -# Atiyah: for more DRY Code, you may consider using -# LPM.ParamTypes.{seed!,reseed0!} within mainHelpers.jl -# and remove the following call & the using statement +include(simPars.analysisFile) + Random.seed!(simPars.seed) # create model object @@ -24,3 +23,16 @@ const logfile = setupLogging(simPars) @time runModel!(model, simPars, pars, logfile) close(logfile) + +if simPars.dumpAgents + open("agents.txt", "w") do f + saveAgents(f, model.pop) + end +end + +if simPars.dumpHouses + open("houses.txt", "w") do f + saveHouses(f, model.houses) + end +end + diff --git a/mainGui.jl b/mainGui.jl index ea34d5a..725c697 100644 --- a/mainGui.jl +++ b/mainGui.jl @@ -1,15 +1,72 @@ -using Raylib -using Raylib: rayvector +include("mainHelpers.jl") -# make this less annoying -const RL = Raylib +include("analysis.jl") -include("lpm.jl") +using GLMakie -include("src/RayGUI/render.jl") -include("src/RayGUI/SimpleGraph.jl") -using .SimpleGraph +struct LTTicks{T} + ticks::T + offset::Float64 + slope::Float64 +end + +function Makie.get_ticks(lt::LTTicks, scale, formatter, vmin, vmax) + lims_transformed = (vmin, vmax) .* lt.slope .+ lt.offset + tickvals_transformed, ticklabels = Makie.get_ticks(lt.ticks, scale, formatter, lims_transformed...) + tickvals_untransformed = (tickvals_transformed .- lt.offset) ./ lt.slope + return tickvals_untransformed, ticklabels +end + + +function setto!(a1::AbstractVector, a2::AbstractVector) + resize!(a1, length(a2)) + a1[:] = a2 +end + +const h_red = colorant"red" +const h_black = colorant"black" + + + +function house_colors!(cols, houses) + resize!(cols, length(houses)) + for (i, h) in enumerate(houses) + cols[i] = isempty(h.occupants) ? h_black : h_red + end +end + +t_pos(t) = t.pos[1]*27, (12-t.pos[2])*27 +h_pos(h) = t_pos(h.town) .+ h.pos + +coords(agent) = h_pos(agent.pos) + +function network!(poss, agent) + empty!.(poss) + ac = coords(agent) + for c in children(agent) + if !alive(c) + continue + end + push!(poss[1], ac) + push!(poss[1], coords(c)) + end + for c in parents(agent) + if isUndefined(c) || !alive(c) + continue + end + push!(poss[2], ac) + push!(poss[2], coords(c)) + end + for c in siblings(agent) + if isUndefined(c) || !alive(c) + continue + end + push!(poss[3], ac) + push!(poss[3], coords(c)) + end +end + function main(parOverrides...) args = copy(ARGS) @@ -19,82 +76,168 @@ function main(parOverrides...) end # need to do that first, otherwise it blocks the GUI - simPars, pars, args = loadParameters(args, - ["--gui-scale"], - Dict(:help => "set gui scale", :default => 1.0, :arg_type => Float64)) + simPars, pars, args = loadParameters(args) + model = setupModel(pars) logfile = setupLogging(simPars) - scale = args[:gui_scale] - screenWidth = floor(Int, 1600 * scale) - screenHeight = floor(Int, 900 * scale) - - RL.InitWindow(screenWidth, screenHeight, "this is a test") - RL.SetTargetFPS(30) - camera = RL.RayCamera2D( - rayvector(screenWidth/2, screenHeight/2), - rayvector(screenWidth/2, screenHeight/2), - #rayvector(500, 500), - 0, - 1 - ) - - # create graph objects with colour - graph_pop = Graph{Float64}(RL.BLUE) - graph_hhs = Graph{Float64}(RL.WHITE) - graph_marr = Graph{Float64}(RL.BLACK) - graph_age = Graph{Float64}(RL.RED) - - pause = false - time = Rational(simPars.startTime) - while !RL.WindowShouldClose() - - if !pause && time <= simPars.finishTime - stepModel!(model, time, simPars, pars) + GLMakie.activate!() + fig = Figure(resolution=(1600,900)) + + #ax_left = Axis(fig[1:2, 1]) + towns = [Rect(t_pos(t)..., 26, 26) for t in model.towns] + cols = [(isempty(t.houses) ? colorant"blue" : colorant"green") for t in model.towns] + ax_left, _ = poly(fig[1:2, 1], towns, color = cols) + hidespines!(ax_left) + hidedecorations!(ax_left) + + houses = [h_pos(h) for h in model.houses] + obs_hc = Observable([h_red]) + house_colors!(obs_hc[], model.houses) + scatter!(houses, color=obs_hc, marker=:rect, markersize=1, markerspace=:data) + + poss = Vector{Vector{Tuple{Int, Int}}}() + push!(poss, []) + push!(poss, []) + push!(poss, []) + obs_poss_c = Observable(poss[1]) + obs_poss_p = Observable(poss[2]) + obs_poss_s = Observable(poss[3]) + lines!(obs_poss_c, color=:yellow) + lines!(obs_poss_p, color=:black) + lines!(obs_poss_s, color=:blue) + f_agent = rand(model.pop) + + dat_pop = [0.0] + dat_marr = [0.0] + obs_pop = Observable([dat_pop, dat_marr]) + ax_pop, _ = series(fig[1,2], obs_pop, labels=["population size", "#married"], + axis=(; xticks=LTTicks(WilkinsonTicks(5), 1920.0, 1/12))) + axislegend(ax_pop) + + obs_age = Observable([0.0]) + ax_age, _ = barplot(fig[1,3], obs_age, label="population pyramid", direction=:x, + axis=(; yticks=LTTicks(WilkinsonTicks(10), 0.0, 3.0))) + axislegend(ax_age) + + obs_careneed = Observable([0.0]) + obs_class = Observable([0.0]) + obs_inc_dec = Observable([0.0]) + obs_age_diff = Observable([0.0]) + ax_careneed, _ = barplot(fig[2,2][1,1], obs_careneed, label="care need") + ax_class, _ = barplot(fig[2,2][1,2], obs_class, label="social class") + ax_inc_dec, _ = barplot(fig[2,2][2,1], obs_inc_dec, label="income deciles") + ax_age_diff, _ = barplot(fig[2,2][2,2], obs_age_diff, label="couple age diff", + axis=(; xticks=LTTicks(WilkinsonTicks(5), -10.0, 1.0))) + axislegend(ax_careneed) + axislegend(ax_class) + axislegend(ax_inc_dec) + axislegend(ax_age_diff) + + dat_cares = [0.0] + dat_careb = [0.0] + obs_care = Observable([dat_cares, dat_careb]) + ax_care, _ = series(fig[2,3], obs_care, labels=["care supply", "unmet care need"], + axis=(; xticks=LTTicks(WilkinsonTicks(5), 1920.0, 1/12))) + axislegend(ax_care) + + #dat_f_status = Vector{Float64}() + #dat_m_status = Vector{Float64}() + + display(fig) + + runbutton = Button(fig[3,3][1,1]; label = "run", tellwidth = false) + pause = Observable(false) + on(runbutton.clicks) do clicks; pause[] = !pause[]; end + quitbutton = Button(fig[3,3][1,2]; label = "quit", tellwidth = false) + goon = Observable(true) + on(quitbutton.clicks) do clicks; goon[]=false; end + + obs_year = Observable("") + Label(fig[3,1][1,1], obs_year, tellwidth=false, fontsize=25) + + randbutton = Button(fig[3,1][1,2]; label = "agent", tellwidth = false) + on(randbutton.clicks) do clicks; f_agent = rand(model.pop); end + + obs_agent= Observable("") + Label(fig[3,1][1,3], obs_agent, tellwidth=false, justification=:left) + + time = Rational(pars.poppars.startTime) + while goon[] + + if !pause[] && time <= pars.poppars.finishTime + stepModel!(model, time, pars) time += simPars.dt - data = observe(Data, model) + data = observe(Data, model, time, pars) log_results(logfile, data) + # add values to graph objects - add_value!(graph_pop, data.alive.n) - add_value!(graph_hhs, data.hh_size.mean) - add_value!(graph_marr, data.married.n) - set_data!(graph_age, data.hist_age.bins, minm=0) - println(data.hh_size.max, " ", data.alive.n, " ", data.eligible.n, " ", data.eligible2.n) + push!(dat_pop, data.alive.n) + push!(dat_marr, data.married.n) + + push!(dat_cares, data.care_supply.mean) + push!(dat_careb, data.unmet_care.mean) + + setto!(obs_careneed[], data.careneed.bins) + setto!(obs_class[], data.class.bins) + setto!(obs_inc_dec[], data.income_deciles) + setto!(obs_age_diff[], data.age_diff.bins) + + setto!(obs_age[], data.age.bins) + + house_colors!(obs_hc[], model.houses) + #setto!(dat_f_status, data.f_status.bins) + #setto!(dat_m_status, data.m_status.bins) + + println(data.hh_size.max, " ", data.alive.n, " ", data.single.n, + " ", data.income.mean) end - - if RL.IsKeyPressed(Raylib.KEY_SPACE) - pause = !pause - sleep(0.2) + + if pause[] + sleep(0.001) end - - RL.BeginDrawing() - - RL.ClearBackground(RL.LIGHTGRAY) + if !alive(f_agent) + f_agent = rand(model.pop) + end + network!(poss, f_agent) - RL.BeginMode2D(camera) + notify(obs_hc) + notify(obs_poss_c) + notify(obs_poss_p) + notify(obs_poss_s) - drawModel(model, (0, 0), - (floor(Int, 50 * scale), floor(Int, 50 * scale)), - (floor(Int, 2 * scale), floor(Int, 2 * scale))) - - RL.EndMode2D() - - # draw graphs - draw_graph(floor(Int, screenWidth/3), 0, floor(Int, screenWidth*2/3), screenHeight, - [graph_pop, graph_hhs, graph_marr, graph_age], - single_scale = false, - labels = ["#alive", "hh size", "#married", "age"], - fontsize = floor(Int, 15 * scale)) + notify(obs_pop) + autolimits!(ax_pop) - - RL.DrawText("$(Float64(time))", 0, - screenHeight - floor(Int, 20 * scale), - floor(Int, 20 * scale), RL.BLACK) - - RL.EndDrawing() + notify(obs_care) + autolimits!(ax_care) + + notify(obs_careneed) + autolimits!(ax_careneed) + notify(obs_class) + autolimits!(ax_class) + notify(obs_inc_dec) + autolimits!(ax_inc_dec) + notify(obs_age_diff) + autolimits!(ax_age_diff) + + notify(obs_age) + autolimits!(ax_age) + + m_status = isUndefined(partner(f_agent)) ? "single" : "married" + m_s = isUndefined(mother(f_agent)) ? "" : "mother" + f_s = isUndefined(father(f_agent)) ? "" : "father" + n_sibs = count(x->!isUndefined(x), siblings(f_agent)) + n_ch = count(x->!isUndefined(x), children(f_agent)) + obs_agent[] = "age: $(floor(Int, age(f_agent)))\n" * + "status: $m_status\n" * + "living parents: $m_s $f_s\n" * + "$n_sibs siblings\n" * + "$n_ch children\n" * + "care need: $(careNeedLevel(f_agent))" + obs_year[] = "$(floor(Int, Float64(time)))" end - RL.CloseWindow() close(logfile) end diff --git a/mainHelpers.jl b/mainHelpers.jl index 1195c68..87024d4 100644 --- a/mainHelpers.jl +++ b/mainHelpers.jl @@ -1,203 +1,37 @@ -addToLoadPath!(".", "src") +include("lib/loadLibsPath.jl") -using ArgParse - -using LPM.ParamTypes +addToLoadPath!(String(@__DIR__) * "/.", + String(@__DIR__) * "/lib", + String(@__DIR__) * "/src") -using XAgents - -using LPM.Demography.Create -using LPM.Demography.Initialize -using LPM.Demography.Simulate +using ArgParse using Utilities +using DemographyPars +using DemographyModelEO -# TODO put into module somewhere? -# Atiyah: Suggestion: as it is related to ParamTypes, it fits there -# or another module for Data (though currently -# not that significant amount of code) -include("src/lpm/demography/demographydata.jl") -include("src/handleParams.jl") - -include("analysis.jl") +include("src/demography/data.jl") -mutable struct Model - towns :: Vector{Town} - houses :: Vector{PersonHouse} - pop :: Vector{Person} - - fertility :: Matrix{Float64} - deathFemale :: Matrix{Float64} - deathMale :: Matrix{Float64} -end - -function createDemography!(pars) - ukTowns = createTowns(pars.mappars) +include("src/handleParams.jl") - ukHouses = Vector{PersonHouse}() - # maybe switch using parameter - #ukPopulation = createPopulation(pars.poppars) - ukPopulation = createPyramidPopulation(pars.poppars) - - # Atiyah: For more DRY code, you may want to consider calling - # loadDemographyData(datapars) +function setupModel(pars) datp = pars.datapars dir = datp.datadir - ukDemoData = loadDemographyData(dir * "/" * datp.fertFName, - dir * "/" * datp.deathFFName, - dir * "/" * datp.deathMFName) - - Model(ukTowns, ukHouses, ukPopulation, - ukDemoData.fertility , ukDemoData.deathFemale, ukDemoData.deathMale) -end - - -function initialConnectH!(houses, towns, pars) - newHouses = initializeHousesInTowns(towns, pars) - append!(houses, newHouses) -end - -function initialConnectP!(pop, houses, pars) - assignCouplesToHouses!(pop, houses) -end - - -function initializeDemography!(model, poppars, workpars, mappars) - initialConnectH!(model.houses, model.towns, mappars) - initialConnectP!(model.pop, model.houses, mappars) - - for person in model.pop - initClass!(person, poppars) - initWork!(person, workpars) - end - - nothing -end - - -# Atiyah: remove this for the primative API simulation function -# alivePeople(model) = Iterators.filter(a->alive(a), model.pop) -# data(model) = model - -function stepModel!(model, time, simPars, pars) - # TODO remove dead people? - - # Atiyah: - doDeaths!(model,time,pars) # a possible unified way - # or the primiative-API - # doDeaths!(alivePeople(model), time, data(model), pars.poppars) - - orphans = Iterators.filter(p->selectAssignGuardian(p), model.pop) - applyTransition!(orphans, assignGuardian!, "adoption", time, model, pars) + demoData = loadDemographyData(dir * "/" * datp.iniAgeFName, + dir * "/" * datp.pre51FertFName, + dir * "/" * datp.fertFName, + dir * "/" * datp.pre51DeathsFName, + dir * "/" * datp.deathFFName, + dir * "/" * datp.deathMFName) - # Atiyah: - babies = doBirths!(model,time,pars) - # babies = doBirths!(alivePeople(model), model.pop), time, model, pars.birthpars) + model = createDemographyModel!(demoData, pars) - selected = Iterators.filter(p->selectAgeTransition(p, pars.workpars), model.pop) - applyTransition!(selected, ageTransition!, "age", time, model, pars.workpars) - - selected = Iterators.filter(p->selectWorkTransition(p, pars.workpars), model.pop) - applyTransition!(selected, workTransition!, "work", time, model, pars.workpars) - - selected = Iterators.filter(p->selectSocialTransition(p, pars.workpars), model.pop) - applyTransition!(selected, socialTransition!, "social", time, model, pars.workpars) - - selected = Iterators.filter(p->selectDivorce(p, pars), model.pop) - applyTransition!(selected, divorce!, "divorce", time, model, - fuse(pars.divorcepars, pars.workpars)) - - resetCacheMarriages() - selected = Iterators.filter(p->selectMarriage(p, pars.workpars), model.pop) - applyTransition!(selected, marriage!, "marriage", time, model, - fuse(pars.poppars, pars.marriagepars, pars.birthpars, pars.mappars)) - - append!(model.pop, babies) -end - - -function loadParameters(argv, cmdl...) - arg_settings = ArgParseSettings("run simulation", autofix_names=true) - - @add_arg_table! arg_settings begin - "--par-file", "-p" - help = "parameter file" - default = "" - "--par-out-file", "-P" - help = "file name for parameter output" - default = "parameters.run.yaml" - end - - if ! isempty(cmdl) - add_arg_table!(arg_settings, cmdl...) - end - - # setup command line arguments with docs + initializeDemographyModel!(model, pars.poppars, pars.workpars, pars.mappars) - add_arg_group!(arg_settings, "Simulation Parameters") - fieldsAsArgs!(arg_settings, SimulationPars) - - for t in fieldtypes(DemographyPars) - groupName = String(nameOfParType(t)) * " Parameters" - add_arg_group!(arg_settings, groupName) - fieldsAsArgs!(arg_settings, t) - end - - # parse command line - args = parse_args(argv, arg_settings, as_symbols=true) - - # read parameters from file if provided or set to default - simpars, pars = loadParametersFromFile(args[:par_file]) - - # override values that were provided on command line - - overrideParsCmdl!(simpars, args) - - @assert typeof(pars) == DemographyPars - for f in fieldnames(DemographyPars) - overrideParsCmdl!(getfield(pars, f), args) - end - - # Atiyah: for more DRY Code, you may consider using - # LPM.ParamTypes.{seed!,reseed0!} within mainHelpers.jl - # and remove the following call & the using statement - # set time dependent seed - if simpars.seed == 0 - simpars.seed = floor(Int, time()) - end - - # keep a record of parameters used (including seed!) - saveParametersToFile(simpars, pars, args[:par_out_file]) - - simpars, pars, args -end - - -function setupModel(pars) - model = createDemography!(pars) - - initializeDemography!(model, pars.poppars, pars.workpars, pars.mappars) - - #= - Atiyah: this portion is not need any more. But if useful, could be - attached with Logging level. - @show "Town Samples: \n" - @show model.towns[1:10] - println(); println(); - - @show "Houses samples: \n" - @show model.houses[1:10] - println(); println(); - - @show "population samples : \n" - @show model.pop[1:10] - println(); println(); - =# - model end @@ -216,20 +50,23 @@ end function runModel!(model, simPars, pars, logfile = nothing; FS = "\t") - time = simPars.startTime + curTime = pars.poppars.startTime simPars.verbose ? setVerbose!() : unsetVerbose!() setDelay!(simPars.sleeptime) - while time < simPars.finishTime - stepModel!(model, time, simPars, pars) + # no point in continuing with the simulation if we are not recording results + finishTime = min(pars.poppars.finishTime, simPars.endLogTime) + + while curTime <= finishTime + stepModel!(model, curTime, pars) - if logfile != nothing - results = observe(Data, model) + if logfile != nothing && curTime >= simPars.startLogTime + results = observe(Data, model, curTime, pars) log_results(logfile, results; FS) end - time += simPars.dt + curTime += simPars.dt end end diff --git a/mainMA.jl b/mainMA.jl deleted file mode 100644 index bb21afe..0000000 --- a/mainMA.jl +++ /dev/null @@ -1,49 +0,0 @@ -""" -Main simulation of the lone parent model - -Run this script from shell as -# julia mainMA.jl - -from REPL execute it using -> include("mainMA.jl") -""" - -include("mainMAHelpers.jl") - -mainConfig = Light() # no input files, logging or flags (REPL Exec.) -# mainConfig = WithInputFiles() - -# lpmExample = LPMUKDemography() # remove deads -lpmExample = LPMUKDemographyOpt() # don't remove deads - -const simPars, pars = loadParameters(mainConfig) - -# Most significant simulation and model parameters -# The following works only with Light() configuration -# useful when executing from REPL -if mainConfig == Light() - simPars.seed = 0; seed!(simPars) - simPars.verbose = false - simPars.checkassumption = false - simPars.sleeptime = 0 - pars.poppars.initialPop = 500 -end - -const model = setupModel(pars) - -const logfile = setupLogging(simPars,mainConfig) - -const demoData = loadDemographyData(pars.datapars) - -const ukDemography = MAModel(model,pars,demoData) - -const lpmDemographySim = - ABMSimulationP{typeof(simPars)}(simPars,setupEnabled = false) - -setup!(lpmDemographySim,lpmExample) - -# Execution -@time run!(ukDemography,lpmDemographySim,lpmExample) - -closeLogfile(logfile,mainConfig) - diff --git a/mainMAHelpers.jl b/mainMAHelpers.jl deleted file mode 100644 index e57c7ab..0000000 --- a/mainMAHelpers.jl +++ /dev/null @@ -1,48 +0,0 @@ -using Random - -include("./loadLibsPath.jl") - -addToLoadPath!("src") -addToLoadPath!("src/multiagents") -addToLoadPath!("../MultiAgents.jl") - -include("mainHelpers.jl") - -using MultiAgents: initMultiAgents, MAVERSION -initMultiAgents() # reset agents counter -@assert MAVERSION == v"0.3.1" # ensure MultiAgents.jl latest update - - -using LPM.ParamTypes: seed! -using MultiAgents: AbstractMABM, ABMSimulationP -using MultiAgents: run! -using MALPM.Demography: MAModel, LPMUKDemography, LPMUKDemographyOpt -using MALPM.Demography.SimSetup: setup! - -""" -How simulations is to be executed: -- with or without input files, arguments and logging -""" - -abstract type MainSim end -struct WithInputFiles <: MainSim end # Input parameter files -struct Light <: MainSim end # no input files - -function loadParameters(::WithInputFiles) - simPars, pars = loadParameters(ARGS) - seed!(simPars) - simPars, pars -end - -function loadParameters(::Light) - simPars = SimulationPars() - seed!(simPars) - pars = DemographyPars() - simPars, pars -end - -setupLogging(simPars,::WithInputFiles) = setupLogging(simPars) -setupLogging(simPars,::Light) = nothing - -closeLogfile(loffile,::WithInputFiles) = close(logfile) -closeLogfile(logfile,::Light) = nothing \ No newline at end of file diff --git a/profile.jl b/profile.jl new file mode 100644 index 0000000..1fc3f18 --- /dev/null +++ b/profile.jl @@ -0,0 +1,29 @@ +using ProfileView + +include("mainHelpers.jl") +include("analysis.jl") + +const simPars, pars = loadParameters(ARGS) + +const model = setupModel(pars) + +t = Rational(pars.poppars.startTime) + +function run_until(t_end) + while true + global model, t, pars, simPars + stepModel!(model, t, pars) + t += simPars.dt + if t >= t_end + break + end + end +end + + +function run_steps(n) + for i in 1:n + global model, t, pars, simPars + stepModel!(model, t, pars) + end +end diff --git a/src/DemographyModel.jl b/src/DemographyModel.jl new file mode 100644 index 0000000..157663e --- /dev/null +++ b/src/DemographyModel.jl @@ -0,0 +1,183 @@ +module DemographyModel + +export Model, createDemographyModel!, initializeDemographyModel!, stepModel! + +include("agents/town.jl") +include("agents/house.jl") +include("agents/person.jl") +include("agents/world.jl") + +include("demography/setup/map.jl") +include("demography/setup/population.jl") +include("demography/setup/mapPop.jl") + +include("demography/simulate/allocate.jl") +include("demography/simulate/death.jl") +include("demography/simulate/birth.jl") +include("demography/simulate/divorce.jl") +include("demography/simulate/ageTransition.jl") +include("demography/simulate/socialTransition.jl") +include("demography/simulate/relocate.jl") +include("demography/simulate/marriages.jl") +include("demography/simulate/dependencies.jl") +include("demography/simulate/socialCareTransition.jl") +include("demography/simulate/care.jl") + +using Utilities + + +mutable struct Model + towns :: Vector{PersonTown} + houses :: Vector{PersonHouse} + pop :: Vector{Person} + babies :: Vector{Person} + + fertFByAge51 :: Vector{Float64} + fertility :: Matrix{Float64} + pre51Fertility :: Vector{Float64} + pre51Deaths :: Matrix{Float64} + deathFemale :: Matrix{Float64} + deathMale :: Matrix{Float64} + + birthCache :: BirthCache{Person} + deathCache :: DeathCache + marriageCache :: MarriageCache{Person} + socialCache :: SocialCache + socialCareCache :: SocialCareCache + divorceCache :: DivorceCache +end + + +function createDemographyModel!(data, pars) + towns = createTowns(pars.mappars) + + houses = Vector{PersonHouse}() + + # maybe switch using parameter + #ukPopulation = createPopulation(pars.poppars) + population = createPyramidPopulation(pars.poppars, data.initialAgePyramid) + + yearsFert = [1951 > data.pre51Fertility[y, 1] >= pars.poppars.startTime + for y in 1:size(data.pre51Fertility)[1]] + + yearsMort = [1951 > data.pre51Deaths[y, 1] >= pars.poppars.startTime + for y in 1:size(data.pre51Deaths)[1]] + + fert = data.fertility[:, 1] # age-specific fertility in 1951 + byAgeF = fert ./ (sum(fert)/length(fert)) + + Model(towns, houses, population, [], + byAgeF, data.fertility, data.pre51Fertility[yearsFert, 2], + data.pre51Deaths[yearsMort, 2:3], data.deathFemale, data.deathMale, + BirthCache{Person}(), DeathCache(), MarriageCache{Person}(), SocialCache(), + SocialCareCache(), DivorceCache()) +end + + +function initialConnectH!(houses, towns, pars) + newHouses = initializeHousesInTowns(towns, pars) + append!(houses, newHouses) +end + +function initialConnectP!(pop, houses, pars) + assignCouplesToHouses!(pop, houses) +end + + +function initializeDemographyModel!(model, poppars, workpars, mappars) + initialConnectH!(model.houses, model.towns, mappars) + initialConnectP!(model.pop, model.houses, mappars) + + for person in model.pop + initClass!(person, poppars) + initWork!(person, workpars) + end + + nothing +end + + +function removeDead!(model) + for i in length(model.pop):-1:1 + if !alive(model.pop[i]) + remove_unsorted!(model.pop, i) + end + end +end + + +function addBaby!(model, baby) + push!(model.babies, baby) +end + + +# TODO not entirely sure if this really belongs here +function stepModel!(model, time, pars) + shuffle!(model.pop) + socialPreCalc!(model, pars) + socialCarePreCalc!(model, fuse(pars.poppars, pars.carepars)) + divorcePreCalc!(model, fuse(pars.poppars, pars.divorcepars, pars.workpars)) + birthPreCalc!(model, fuse(pars.poppars, pars.birthpars)) + deathPreCalc!(model, pars.poppars) + + applyTransition!(model.pop, "death") do person + death!(person, time, model, pars.poppars) + end + removeDead!(model) + + orphans = Iterators.filter(p->selectAssignGuardian(p), model.pop) + applyTransition!(orphans, "adoption") do person + assignGuardian!(person, time, model, pars) + end + + selected = Iterators.filter(p->selectBirth(p, pars.birthpars), model.pop) + applyTransition!(selected, "birth") do person + birth!(person, time, model, fuse(pars.poppars, pars.birthpars), addBaby!) + end + + selected = Iterators.filter(p->selectAgeTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "age") do person + ageTransition!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectSocialCareTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "social care") do person + socialCareTransition!(person, time, model, fuse(pars.poppars, pars.carepars)) + end + + socialCare!(model, pars.carepars) + + selected = Iterators.filter(p->selectWorkTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "work") do person + workTransition!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectRelocate(p, pars.workpars), model.pop) + applyTransition!(selected, "relocate") do person + relocate!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectSocialTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "social") do person + socialTransition!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectDivorce(p, pars), model.pop) + applyTransition!(selected, "divorce") do person + divorce!(person, time, model, fuse(pars.poppars, pars.divorcepars, pars.workpars)) + end + + marriagePreCalc!(model, fuse(pars.poppars, pars.marriagepars, pars.birthpars, pars.mappars)) + selected = Iterators.filter(p->selectMarriage(p, pars.workpars), model.pop) + applyTransition!(selected, "marriage") do person + marriage!(person, time, model, + fuse(pars.poppars, pars.marriagepars, pars.birthpars, pars.mappars)) + end + + + append!(model.pop, model.babies) + empty!(model.babies) +end + + +end # module DemographyModel diff --git a/src/DemographyModelEO.jl b/src/DemographyModelEO.jl new file mode 100644 index 0000000..ee1f388 --- /dev/null +++ b/src/DemographyModelEO.jl @@ -0,0 +1,183 @@ +module DemographyModelEO + +export Model, createDemographyModel!, initializeDemographyModel!, stepModel! + +include("agents/town.jl") +include("agents/house.jl") +include("agents/person.jl") +include("agents/world.jl") + +include("demography/setup/map.jl") +include("demography/setup/population.jl") +include("demography/setup/mapPop.jl") + +include("demography/simulate/allocate.jl") +include("demography/simulate/death.jl") +include("demography/simulate/birth.jl") +include("demography/simulate/divorce.jl") +include("demography/simulate/ageTransition.jl") +include("demography/simulate/socialTransition.jl") +include("demography/simulate/relocate.jl") +include("demography/simulate/marriages.jl") +include("demography/simulate/dependencies.jl") +include("demography/simulate/socialCareTransition.jl") +include("demography/simulate/care.jl") + +using Utilities + + +mutable struct Model + towns :: Vector{PersonTown} + houses :: Vector{PersonHouse} + pop :: Vector{Person} + babies :: Vector{Person} + + fertFByAge51 :: Vector{Float64} + fertility :: Matrix{Float64} + pre51Fertility :: Vector{Float64} + pre51Deaths :: Matrix{Float64} + deathFemale :: Matrix{Float64} + deathMale :: Matrix{Float64} + + birthCache :: BirthCache{Person} + deathCache :: DeathCache + marriageCache :: MarriageCache{Person} + socialCache :: SocialCache + socialCareCache :: SocialCareCache + divorceCache :: DivorceCache +end + + +function createDemographyModel!(data, pars) + towns = createTowns(pars.mappars) + + houses = Vector{PersonHouse}() + + # maybe switch using parameter + #ukPopulation = createPopulation(pars.poppars) + population = createPyramidPopulation(pars.poppars, data.initialAgePyramid) + + yearsFert = [1951 > data.pre51Fertility[y, 1] >= pars.poppars.startTime + for y in 1:size(data.pre51Fertility)[1]] + + yearsMort = [1951 > data.pre51Deaths[y, 1] >= pars.poppars.startTime + for y in 1:size(data.pre51Deaths)[1]] + + fert = data.fertility[:, 1] # age-specific fertility in 1951 + byAgeF = fert ./ (sum(fert)/length(fert)) + + Model(towns, houses, population, [], + byAgeF, data.fertility, data.pre51Fertility[yearsFert, 2], + data.pre51Deaths[yearsMort, 2:3], data.deathFemale, data.deathMale, + BirthCache{Person}(), DeathCache(), MarriageCache{Person}(), SocialCache(), + SocialCareCache(), DivorceCache()) +end + + +function initialConnectH!(houses, towns, pars) + newHouses = initializeHousesInTowns(towns, pars) + append!(houses, newHouses) +end + +function initialConnectP!(pop, houses, pars) + assignCouplesToHouses!(pop, houses) +end + + +function initializeDemographyModel!(model, poppars, workpars, mappars) + initialConnectH!(model.houses, model.towns, mappars) + initialConnectP!(model.pop, model.houses, mappars) + + for person in model.pop + initClass!(person, poppars) + initWork!(person, workpars) + end + + nothing +end + + +function removeDead!(model) + for i in length(model.pop):-1:1 + if !alive(model.pop[i]) + remove_unsorted!(model.pop, i) + end + end +end + + +function addBaby!(model, baby) + push!(model.babies, baby) +end + + +# TODO not entirely sure if this really belongs here +function stepModel!(model, time, pars) + shuffle!(model.pop) + socialPreCalc!(model, pars) + socialCarePreCalc!(model, fuse(pars.poppars, pars.carepars)) + divorcePreCalc!(model, fuse(pars.poppars, pars.divorcepars, pars.workpars)) + birthPreCalc!(model, fuse(pars.poppars, pars.birthpars)) + deathPreCalc!(model, pars.poppars) + + selected = Iterators.filter(p->selectAgeTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "age") do person + ageTransition!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectBirth(p, pars.birthpars), model.pop) + applyTransition!(selected, "birth") do person + birth!(person, time, model, fuse(pars.poppars, pars.birthpars), addBaby!) + end + + orphans = Iterators.filter(p->selectAssignGuardian(p), model.pop) + applyTransition!(orphans, "adoption") do person + assignGuardian!(person, time, model, pars) + end + + selected = Iterators.filter(p->selectWorkTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "work") do person + workTransition!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectSocialTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "social") do person + socialTransition!(person, time, model, pars.workpars) + end + + marriagePreCalc!(model, fuse(pars.poppars, pars.marriagepars, pars.birthpars, pars.mappars)) + selected = Iterators.filter(p->selectMarriage(p, pars.workpars), model.pop) + applyTransition!(selected, "marriage") do person + marriage!(person, time, model, + fuse(pars.poppars, pars.marriagepars, pars.birthpars, pars.mappars)) + end + + selected = Iterators.filter(p->selectRelocate(p, pars.workpars), model.pop) + applyTransition!(selected, "relocate") do person + relocate!(person, time, model, pars.workpars) + end + + selected = Iterators.filter(p->selectDivorce(p, pars), model.pop) + applyTransition!(selected, "divorce") do person + divorce!(person, time, model, fuse(pars.poppars, pars.divorcepars, pars.workpars)) + end + + selected = Iterators.filter(p->selectSocialCareTransition(p, pars.workpars), model.pop) + applyTransition!(selected, "social care") do person + socialCareTransition!(person, time, model, fuse(pars.poppars, pars.carepars)) + end + + socialCare!(model, pars.carepars) + + applyTransition!(model.pop, "death") do person + death!(person, time, model, pars.poppars) + end + removeDead!(model) + + + append!(model.pop, model.babies) + empty!(model.babies) +end + + +end # module DemographyModel diff --git a/src/lpm/ParamTypes.jl b/src/DemographyPars.jl similarity index 76% rename from src/lpm/ParamTypes.jl rename to src/DemographyPars.jl index ab64fc3..0705748 100644 --- a/src/lpm/ParamTypes.jl +++ b/src/DemographyPars.jl @@ -1,24 +1,20 @@ """ Parameter types and values - -This module is within the LPM module """ -module ParamTypes +module DemographyPars using Random using Parameters import Random.seed! -export SimulationPars, seed!, reseed0! +export SimulationPars, reseed0!, seed! -include("./demography/demographypars.jl") +include("./demography/parameters.jl") "General simulation parameters" @with_kw mutable struct SimulationPars - startTime :: Rational{Int} = 1920 - finishTime :: Rational{Int} = 2040 dt :: Rational{Int} = 1//12 # step size seed :: Int = 42 ; @assert seed >= 0 # 0 is random verbose::Bool = false # whether significant intermediate info shallo be printed @@ -26,6 +22,12 @@ include("./demography/demographypars.jl") # how long simulation is suspended after printing info checkassumption :: Bool = false # whether assumptions in unit functions should be checked logfile :: String = "log.tsv" + analysisFile :: String = "analysis.jl" + startLogTime :: Rational{Int} = 0 + endLogTime :: Rational{Int} = 10000 + + dumpAgents :: Bool = false + dumpHouses :: Bool = false end reseed0!(simPars) = @@ -41,4 +43,4 @@ function seed!(simPars::SimulationPars, end -end # ParamTypes +end # DemographyPars diff --git a/src/LPM.jl b/src/LPM.jl deleted file mode 100644 index b376505..0000000 --- a/src/LPM.jl +++ /dev/null @@ -1,9 +0,0 @@ -module LPM - export LPMVERSION - - const LPMVERSION = r"0.6.1" - - include("./lpm/ParamTypes.jl") - include("./lpm/Demography.jl") - -end # LoneParentsModel diff --git a/src/MALPM.jl b/src/MALPM.jl deleted file mode 100644 index 83ccc11..0000000 --- a/src/MALPM.jl +++ /dev/null @@ -1,8 +0,0 @@ -""" -Main components of LoneParentsModel simulation based on MultiAgents.jl package -""" -module MALPM - - include("./malpm/Demography.jl") - -end # MALPM \ No newline at end of file diff --git a/src/RayGUI/SimpleGraph.jl b/src/RayGUI/SimpleGraph.jl deleted file mode 100644 index 4a211a9..0000000 --- a/src/RayGUI/SimpleGraph.jl +++ /dev/null @@ -1,123 +0,0 @@ -# Copyright (C) 2020 Martin Hinsch -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - - - -module SimpleGraph - -using Printf - -using Raylib -const RL = Raylib - -export Graph, add_value!, draw_graph, set_data! - -### super simplistic graph implementation - -mutable struct Graph{T} - data :: Vector{T} - max :: T - min :: T - colour :: RL.RayColor -end - -Graph{T}(col) where {T} = Graph{T}([], typemin(T), typemax(T), col) - -function add_value!(graph::Graph, value) - push!(graph.data, value) - value > graph.max ? (graph.max = value) : (value < graph.min ? (graph.min = value) : value) -end - -function set_data!(graph::Graph, data; maxm = data[1], minm = data[1]) - graph.data = data - graph.max = maxm == data[1] ? maximum(data) : maxm - graph.min = minm == data[1] ? minimum(data) : minm -end - - -function draw_legend(value, x, y, fontsize, colour) - text = @sprintf "%g" value - RL.DrawText(text, x, y, fontsize, colour) -end - - -# draw graph to canvas -function draw_graph(x_0, y_0, width, height, graphs; - single_scale=true, labels=[], fontsize=15) - if single_scale # draw all graphs to the same y scale - max_all = mapreduce(g -> g.max, max, graphs) # find maximum of graphs[...].max - min_all = mapreduce(g -> g.min, min, graphs) - end - - width_legend = RL.MeasureText("00000000", fontsize) - width_g = width - width_legend - height_g = height - fontsize - x_0_g = x_0 + width_legend - - for g in graphs - g_max = single_scale ? max_all : g.max - g_min = single_scale ? min_all : g.min - - # no x or y range, can't draw - if g_max <= g_min || length(g.data) <= 1 - continue - end - - x_scale = (width_g-1) / (length(g.data)-1) - y_scale = (height_g-1) / (g_max - g_min) - - dxold = 1 - dyold = height_g - trunc(Int, (g.data[1]-g_min) * y_scale ) - - for (i,dat) in enumerate(g.data) - dx = trunc(Int, (i-1) * x_scale) + 1 - dy = height_g - trunc(Int, (dat-g_min) * y_scale) - RL.DrawLine(x_0_g+dxold, y_0+dyold, x_0_g+dx, y_0+dy, g.colour) - dxold, dyold = dx, dy - end - end - - if single_scale - draw_legend(min_all, x_0, y_0 + height, fontsize, graphs[1].colour) - draw_legend(max_all, x_0, y_0, fontsize, graphs[1].colour) - else - yoffs = y_0 + height - fontsize * length(graphs) - for (i, g) in enumerate(graphs) - draw_legend(g.min, x_0, yoffs + (i-1) * fontsize, fontsize, g.colour) - draw_legend(g.max, x_0, y_0 + (i-1) * fontsize, fontsize, g.colour) - end - end - - if isempty(labels) - return nothing - end - - @assert length(labels) == length(graphs) - - w = 0 - for l in labels - w = max(w, RL.MeasureText(l*" ", fontsize)) - end - - lx = x_0 + width - w - - for (i, l) in enumerate(labels) - RL.DrawText(l, lx, (i-1)*fontsize, fontsize, graphs[i].colour) - end - -end - - -end diff --git a/src/RayGUI/render.jl b/src/RayGUI/render.jl deleted file mode 100644 index b9f404f..0000000 --- a/src/RayGUI/render.jl +++ /dev/null @@ -1,25 +0,0 @@ - - -function drawModel(model, offset, townSize, houseSize) - for house in model.houses - town_pos = house.town.pos - # top left corner of town - t_offset = offset[1] + (town_pos[1]-1) * townSize[1], - offset[2] + (town_pos[2]-1) * townSize[2] - - # top left corner of house - h_offset = t_offset[1] + (house.pos[1]-1) * houseSize[1], - t_offset[2] + (house.pos[2]-1) * houseSize[2] - - col = isempty(house.occupants) ? RL.BLACK : RL.RED - RL.DrawRectangleRec(RL.RayRectangle(h_offset..., houseSize...), col) - end - - for town in model.towns - t_offset = offset[1] + (town.pos[1]-1) * townSize[1], - offset[2] + (town.pos[2]-1) * townSize[2] - - RL.DrawRectangleLinesEx(RL.RayRectangle(t_offset..., townSize...), 1.0, RL.GREEN) - end -end - diff --git a/src/Utilities.jl b/src/Utilities.jl index ba9a724..e117b19 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -3,42 +3,23 @@ Diverse useful functions and types """ module Utilities -# Types -export Gender, male, female, unknown -# Constants -export SimulationFolderPrefix # Functions -export createTimeStampedFolder, p_yearly2monthly, applyTransition!, remove_unsorted! -export removefirst!, date2yearsmonths, age2yearsmonths +export p_yearly2monthly, applyTransition!, remove_unsorted!, limit +export removefirst!, sorted_unique!, date2yearsmonths, age2yearsmonths export checkAssumptions!, ignoreAssumptions!, assumption, setDelay!, delay export setVerbose!, unsetVerbose!, verbose, verbosePrint, delayedVerbose -export fuse +export fuse, countSubset, @static_var +export dump, dump_property, dump_header +export sumClassBias, rateBias, preCalcRateBias! # list of types -"Gender type enumeration" -@enum Gender unknown female male # constants -"Folder in which simulation results are stored" -const SimulationFolderPrefix = "Simulations_Folder" - -# timeStamp ... - -"create a folder in which simulation results are stored" -function createTimeStampedFolder() - #timeStamp = datetime.datetime.today().strftime('%Y_%m_%d-%H_%M_%S') - #folder = os.path.join('Simulations_Folder', timeStamp) - #if not os.path.exists(folder): - # os.makedirs(folder) - # folder - "" -end - "remove first occurance of e in list" function removefirst!(list, e) e ∉ list ? throw(ArgumentError("element $(e) not in $(list)")) : nothing @@ -59,22 +40,7 @@ age2yearsmonths(age) = date2yearsmonths(age) p_yearly2monthly(p) = 1 - (1-p)^(1/12) -# constants - -"Folder in which simulation results are stored" -const SimulationFolderPrefix = "Simulations_Folder" - -# timeStamp ... - -"create a folder in which simulation results are stored" -function createTimeStampedFolder() - #timeStamp = datetime.datetime.today().strftime('%Y_%m_%d-%H_%M_%S') - #folder = os.path.join('Simulations_Folder', timeStamp) - #if not os.path.exists(folder): - # os.makedirs(folder) - # folder - "" -end +limit(mi, v, ma) = min(ma, max(mi, v)) "Very efficiently remove element `index` from `list`. Does not preserve ordering of `list`." function remove_unsorted!(list, index) @@ -82,11 +48,25 @@ function remove_unsorted!(list, index) pop!(list) end +function sorted_unique!(ar) + v = ar[end] + l = max(0, length(ar) - 1) + for i in l:-1:1 + e = ar[i] + if e == v + ar[i] = ar[end] + pop!(ar) + else + v = e + end + end +end + "Apply a transition function to an iterator." -function applyTransition!(people, transition, name, args...) +function applyTransition!(transition, people, name) count = 0 for p in people - transition(p, args...) + transition(p) count += 1 end @@ -97,7 +77,53 @@ function applyTransition!(people, transition, name, args...) end end +"keep variable across function calls" +macro static_var(init) + var = gensym() + Base.eval(__module__, :(const $var = $init)) + quote + global $var + $var + end |> esc +end + + +"Count the elements of a subset and a subset of that subset of population." +@inline function countSubset(condAll, condSubset, population) + nAll = 0 + nS = 0 + + for x in Iterators.filter(condAll, population) + nAll += 1 + if condSubset(x) + nS += 1 + end + end + + nAll, nS +end + +@inline function sumClassBias(classFn, range, bias) + sum(range) do c + classFn(c) * bias^c + end +end + +@inline function rateBias(classFn, range, bias, thisClass) + bias^thisClass / sumClassBias(classFn, range, bias) +end + + +@inline function preCalcRateBias!(fn, range, bias, array, offset = 0) + sumBias = sumClassBias(fn, range, bias) + for c in range + array[c+offset] = bias^c/sumBias + end +end + + + mutable struct Debug checkAssumptions :: Bool verbose :: Bool @@ -161,4 +187,59 @@ end # both put together :($tuptyp($tup)) end + + +function dump_header(io, obj, FS="\t") + fns = join(fieldnames(typeof(obj)), FS) + print(io, fns) +end + +function dump_property(io, prop, FS="\t", ES=",") + print(io, prop) +end + +function dump_property(io, prop::Rational, FS="\t", ES=",") + print(io, Float64(prop)) +end + +function dump_property(io, prop::Vector, FS="\t", ES=",") + print(io, "(") + for (i, el) in enumerate(prop) + dump_property(io, el, FS, ES) + if i != length(prop) + print(io, ES) + end + end + print(io, ")") +end + +function dump_property(io, prop::Matrix, FS="\t", ES=",") + print(io, "(") + for i in 1:size(prop)[1] + print(io, "(") + for j in 1:size(prop)[2] + el = prop[i, j] + dump_property(io, el, FS, ES) + if j != size(prop)[2] + print(io, ES) + end + end + print(io, ")") + if i != size(prop)[1] + print(io, ES) + end + end + print(io, ")") +end + +function dump(io, obj, FS="\t", ES=",") + for (i, f) in enumerate(fieldnames(typeof(obj))) + dump_property(io, getfield(obj, f), FS, ES) + if i != fieldcount(typeof(obj)) + print(io, FS) + end + end +end + + end # module Utilities diff --git a/src/agents/agents_modules/basicinfo.jl b/src/agents/agents_modules/basicinfo.jl index 560a681..f84076b 100644 --- a/src/agents/agents_modules/basicinfo.jl +++ b/src/agents/agents_modules/basicinfo.jl @@ -1,43 +1,29 @@ -using Utilities: Gender, unknown, female, male, age2yearsmonths +using Utilities export isFemale, isMale, agestep!, agestepAlive!, hasBirthday, yearsold +export Gender, male, female, unknown + +"Gender type enumeration" +@enum Gender unknown female male # TODO think about whether to make this immutable -mutable struct BasicInfoBlock - age::Rational{Int} +@kwdef struct BasicInfo + alive :: Bool = true + age::Rational{Int} = 0//1 # (birthyear, birthmonth) - gender::Gender - alive::Bool -end - -"Default constructor" -BasicInfoBlock(;age=0//1, gender=unknown, alive = true) = BasicInfoBlock(age,gender,alive) - -isFemale(person::BasicInfoBlock) = person.gender == female -isMale(person::BasicInfoBlock) = person.gender == male - - -"costum @show method for Agent person" -function Base.show(io::IO, info::BasicInfoBlock) - year, month = age2yearsmonths(info.age) - print(" $(year) years & $(month) months, $(info.gender) ") - info.alive ? print(" alive ") : print(" dead ") + gender::Gender = unknown end -# setDead!(person::BasicInfoBlock) = person.alive = false +isFemale(person) = person.gender == female +isMale(person) = person.gender == male "increment an age for a person to be used in typical stepping functions" -agestep!(person::BasicInfoBlock, dt=1//12) = person.age += dt +agestep!(person, dt=1//12) = person.age += dt -"increment an age for a person to be used in typical stepping functions" -function agestepAlive!(person::BasicInfoBlock, dt=1//12) - person.age += person.alive ? dt : 0 -end +hasBirthday(person) = person.age % 1 == 0 -hasBirthday(person::BasicInfoBlock) = person.age % 1 == 0 - -function yearsold(person::BasicInfoBlock) - years, = age2yearsmonths(person.age) - years +function yearsold(person) + trunc(Int, person.age) end + diff --git a/src/agents/agents_modules/care.jl b/src/agents/agents_modules/care.jl index 19b692c..cbb5d32 100644 --- a/src/agents/agents_modules/care.jl +++ b/src/agents/agents_modules/care.jl @@ -1,7 +1,6 @@ -export CareBlock -mutable struct CareBlock - careNeedLevel :: Int - socialWork :: Int - childWork :: Int +@kwdef struct Care + careNeedLevel :: Int = 0 + socialWork :: Int = 0 + childWork :: Int = 0 end diff --git a/src/agents/agents_modules/carehouse.jl b/src/agents/agents_modules/carehouse.jl new file mode 100644 index 0000000..0479266 --- /dev/null +++ b/src/agents/agents_modules/carehouse.jl @@ -0,0 +1,38 @@ +export provideCare!, receiveCare!, resetCare!, careBalance + + +struct CareHouse{H} + "net care this house produces (or demands for values < 0)" + netCareSupply :: Int + "net care this house exports to others (or receives for values < 0)" + careProvided :: Int + careConnections :: Vector{H} +end + +function provideCare!(house, care, to) + house.careProvided += care + if !(to in house.careConnections) + push!(house.careConnections, to) + end + #sort!(house.careConnections, by=objectid) + #sorted_unique!(house.careConnections) + #unique!(house.careConnections) +end + +function receiveCare!(house, care, from) + house.careProvided -= care + if !(from in house.careConnections) + push!(house.careConnections, from) + end + #sort!(house.careConnections, by=objectid) + #sorted_unique!(house.careConnections) + #unique!(house.careConnections) +end + +function resetCare!(house, ncs) + house.netCareSupply = ncs + house.careProvided = 0 + empty!(house.careConnections) +end + +careBalance(house) = house.netCareSupply - house.careProvided diff --git a/src/agents/agents_modules/class.jl b/src/agents/agents_modules/class.jl index 96b8675..169c305 100644 --- a/src/agents/agents_modules/class.jl +++ b/src/agents/agents_modules/class.jl @@ -1,9 +1,8 @@ -export ClassBlock - export addClassRank! -mutable struct ClassBlock - classRank :: Int +@kwdef struct Class + classRank :: Int = 0 + parentClassRank :: Int = 0 end -addClassRank!(class, n) = class.classRank += 1 +addClassRank!(class, n=1) = class.classRank += n diff --git a/src/agents/agents_modules/dependencies.jl b/src/agents/agents_modules/dependencies.jl index dd4fe4e..ccf5eeb 100644 --- a/src/agents/agents_modules/dependencies.jl +++ b/src/agents/agents_modules/dependencies.jl @@ -1,15 +1,10 @@ -export DependencyBlock - - -mutable struct DependencyBlock{P} - guardians :: Vector{P} - dependents :: Vector{P} - provider :: Union{P, Nothing} - providees :: Vector{P} +@kwdef struct Dependency{P} + guardians :: Vector{P} = [] + dependents :: Vector{P} = [] + provider :: P = undefinedPerson#undefined(P) + providees :: Vector{P} = [] end -DependencyBlock{P}() where {P} = DependencyBlock{P}([], [], nothing, []) - isDependent(p) = !isempty(p.guardians) hasDependents(p) = isempty(p.dependents) diff --git a/src/agents/agents_modules/kinship.jl b/src/agents/agents_modules/kinship.jl index ffdbade..0a42847 100644 --- a/src/agents/agents_modules/kinship.jl +++ b/src/agents/agents_modules/kinship.jl @@ -1,26 +1,32 @@ -export KinshipBlock -export hasChildren, addChild!, isSingle, parents, siblings - -mutable struct KinshipBlock{P} - father::Union{P,Nothing} - mother::Union{P,Nothing} - partner::Union{P,Nothing} - children::Vector{P} +using Utilities + +export hasChildren, addChild!, isSingle, parents, siblings, nChildren + +@kwdef struct Kinship{P} + father::P = undefinedPerson#undefined(P) + mother::P = undefinedPerson#undefined(P) + partner::P = undefinedPerson#undefined(P) + pTime :: Rational{Int} = 0//1 + children::Vector{P} = [] end -hasChildren(parent::KinshipBlock{P}) where{P} = length(parent.children) > 0 +hasChildren(parent) = length(parent.children) > 0 + +function addChild!(parent, child) + push!(parent.children, child) +end -addChild!(parent::KinshipBlock{P}, child::P) where{P} = push!(parent.children, child) +isSingle(person) = isUndefined(person.partner) -isSingle(person::KinshipBlock) = person.partner == nothing +parents(person) = [person.father, person.mother] -parents(person::KinshipBlock) = [person.father, person.mother] +nChildren(person) = length(person.children) -function siblings(person::KinshipBlock{P}) where P +function siblings(person::P) where {P} sibs = P[] for p in parents(person) - if p == nothing continue end + if isUndefined(p) continue end for c in children(p) if c != person push!(sibs, c) @@ -31,16 +37,3 @@ function siblings(person::KinshipBlock{P}) where P sibs end -"costum @show method for Agent person" -function Base.show(io::IO, kinship::KinshipBlock) - father = kinship.father; mother = kinship.mother; partner = kinship.partner; children = kinship.children; - father == nothing ? nothing : print(" , father : $(father.id)") - mother == nothing ? nothing : print(" , mother : $(mother.id)") - partner == nothing ? nothing : print(" , partner : $(partner.id)") - length(children) == 0 ? nothing : print(" , children : ") - for child in children - print(" $(child.id) ") - end - println() -end - diff --git a/src/agents/agents_modules/maternity.jl b/src/agents/agents_modules/maternity.jl index 134c4e5..3573876 100644 --- a/src/agents/agents_modules/maternity.jl +++ b/src/agents/agents_modules/maternity.jl @@ -1,9 +1,8 @@ -export MaternityBlock export startMaternity!, stepMaternity!, endMaternity!, isInMaternity, maternityDuration -mutable struct MaternityBlock - maternityStatus :: Bool - monthsSinceBirth :: Int +@kwdef struct Maternity + maternityStatus :: Bool = false + monthsSinceBirth :: Int = 0 end isInMaternity(mat) = mat.maternityStatus diff --git a/src/agents/agents_modules/work.jl b/src/agents/agents_modules/work.jl index b8b87bd..8288fd6 100644 --- a/src/agents/agents_modules/work.jl +++ b/src/agents/agents_modules/work.jl @@ -1,6 +1,5 @@ using EnumX -export WorkBlock export setEmptyJobSchedule!, setFullWeeklyTime! export WorkStatus @@ -8,36 +7,36 @@ export WorkStatus # better (scoped) enums from package EnumX @enumx WorkStatus child teenager student worker retired unemployed +const WST = WorkStatus.T + # TODO some of this should probably be moved to Care -mutable struct WorkBlock +@kwdef struct Work "work status" - status :: WorkStatus.T - outOfTownStudent :: Bool - newEntrant :: Bool - initialIncome :: Float64 - finalIncome :: Float64 - wage :: Float64 - income :: Float64 + status :: WST = WorkStatus.child + outOfTownStudent :: Bool = false + newEntrant :: Bool = true + initialIncome :: Float64 = 0 + finalIncome :: Float64 = 0 + wage :: Float64 = 0 + income :: Float64 = 0 "income for full work schedule" - potentialIncome :: Float64 + potentialIncome :: Float64 = 0 "periods worked so far in current job" - jobTenure :: Int + jobTenure :: Int = 0 "7x24 schedule of actual working hours" - schedule :: Matrix{Int} + schedule :: Matrix{Int} = zeros(Int, 7, 24) "potential total working hours per week" - workingHours :: Int + workingHours :: Int = 0 "free time slots" - weeklyTime :: Matrix{Int} + weeklyTime :: Matrix{Int} = zeros(Int, 7, 24) "sum of realised working hours" - availableWorkingHours :: Int + availableWorkingHours :: Int = 0 "lifetime work" - workingPeriods :: Int - pension :: Float64 + workingPeriods :: Float64 = 0 + workExperience :: Float64 = 0 + pension :: Float64 = 0 end -WorkBlock() = WorkBlock(WorkStatus.child, false, true, 0, 0, 0, 0, 0, 0, zeros(Int, 7, 24), - 0, zeros(Int, 7, 24), 0, 0, 0) - function setEmptyJobSchedule!(work) work.schedule = zeros(Int, 7, 24) end @@ -45,3 +44,4 @@ end function setFullWeeklyTime!(work) work.weeklyTime = ones(Int, 7, 24) end + diff --git a/src/agents/house.jl b/src/agents/house.jl index 96ba5ec..1a392b8 100644 --- a/src/agents/house.jl +++ b/src/agents/house.jl @@ -1,32 +1,40 @@ +using CompositeStructs + + export House, HouseLocation export getHomeTown, getHouseLocation, undefined, isEmpty, town -using Utilities: removefirst! +using Utilities +using DeclUtils +include("agents_modules/carehouse.jl") -const HouseLocation = NTuple{2,Int} - -""" -Specification of a House Agent Type. -This file is included in the module XAgents +const HouseLocation = NTuple{2,Int} -Type House to extend from AbstracXAgent. -""" -mutable struct House{P, T} <: AbstractXAgent - id :: Int +mutable struct House{P, T} town :: T pos :: HouseLocation # location in the town # size::String # TODO enumeration type / at the moment not yet necessary occupants::Vector{P} - - House{P, T}(town, pos) where {P, T} = new(getIDCOUNTER(),town, pos,P[]) + + # manual composition for now, @composite does not allow + # partial specialisation + + "net care this house produces (or demands for values < 0)" + netCareSupply :: Int + "net care this house exports to others (or receives for values < 0)" + careProvided :: Int + careConnections :: Vector{House{P, T}} end # House +House{P, T}(t, p) where{P, T} = House(t, p, P[], 0, 0, House{P, T}[]) + + +@export_forward House [netCareSupply, careProvided, careConnections] -undefined(house) = house.town == undefinedTown && house.pos == (-1,-1) isEmpty(house) = length(house.occupants) == 0 @@ -34,28 +42,28 @@ town(house) = house.town # to replace the functions below in order to unify style across agents APIs "town associated with house" -function getHomeTown(house::House) +function getHomeTown(house) house.town end "town name associated with house" -function getHomeTownName(house::House) +function getHomeTownName(house) house.town.name end "house location in the associated town" -function getHouseLocation(house::House) +function getHouseLocation(house) house.pos end "add an occupant to a house" -function addOccupant!(house::House{P}, person::P) where {P} +function addOccupant!(house, person) push!(house.occupants, person) nothing end "remove an occupant from a house" -function removeOccupant!(house::House{P}, person::P) where {P} +function removeOccupant!(house, person) removefirst!(house.occupants, person) # we can't assume anything about the layout of typeof(person) #person.pos = undefinedHouse @@ -63,7 +71,7 @@ function removeOccupant!(house::House{P}, person::P) where {P} end "Costum print function for agents" -function Base.show(io::IO, house::House{P}) where P +function Base.show(io::IO, house::House) townName = getHomeTownName(house) print("House $(house.id) @ town $(townName) @ $(house.pos)") length(house.occupants) == 0 ? nothing : print(" occupants: ") @@ -72,3 +80,16 @@ function Base.show(io::IO, house::House{P}) where P end println() end + + +function Utilities.dump_header(io, h::House, FS) + print(io, "id", FS, "pos", FS) + Utilities.dump_header(io, h.care, FS); print(io, FS) +end + +function Utilities.dump(io, house::House, FS="\t", ES=",") + print(io, objectid(house), FS) + Utilities.dump_property(io, house.pos, FS, ES); print(io, FS) + # no need to dump inhabitants as well, they link back anyway + Utilities.dump(io, house.care, FS, ES); print(io, FS) +end diff --git a/src/agents/household.jl b/src/agents/household.jl deleted file mode 100644 index bf0f202..0000000 --- a/src/agents/household.jl +++ /dev/null @@ -1,15 +0,0 @@ -export Household - -""" -Household concept. -It is still not determined if a household is an agent or an ABM -within a multi-ABM. -""" -mutable struct Household <: AbstractXAgent - id # list of person ids - pos # House - income - # ... -end - -# Constructor Household(personList) diff --git a/src/agents/person.jl b/src/agents/person.jl index 9baa6bc..154072d 100644 --- a/src/agents/person.jl +++ b/src/agents/person.jl @@ -1,13 +1,10 @@ -using TypedDelegation +using CompositeStructs using DeclUtils - -include("../../loadLibsPath.jl") -# enable using/import from local directory -addToLoadPath!("agents_modules") +using TypedDelegation export Person -export PersonHouse, undefinedHouse +export PersonHouse, isUndefined, undefinedHouse export moveToHouse!, resetHouse!, resolvePartnership!, householdIncome export householdIncomePerCapita @@ -34,127 +31,64 @@ include("agents_modules/dependencies.jl") """ Specification of a Person Agent Type. - -This file is included in the module XAgents - -Type Person extends from AbstractAgent. - Person ties various agent modules into one compound agent type. """ -# vvv More classification of attributes (Basic, Demography, Relatives, Economy ) -mutable struct Person <: AbstractXAgent - id::Int - """ - location of a parson's house in a map which implicitly - - (x-y coordinates of a house) - - (town::Town, x-y location in the map) - """ - pos::House{Person} - info::BasicInfoBlock - kinship::KinshipBlock{Person} - maternity :: MaternityBlock - work :: WorkBlock - care :: CareBlock - class :: ClassBlock - dependencies :: DependencyBlock{Person} - - # Person(id,pos,age) = new(id,pos,age) - "Internal constructor" - function Person(pos, info, kinship, maternity, work, care, class, dependencies) - person = new(getIDCOUNTER(),pos,info,kinship, maternity, work, care, class, dependencies) - if !undefined(pos) - addOccupant!(pos, person) - end - if kinship.father != nothing - addChild!(kinship.father,person) - end - if kinship.mother != nothing - addChild!(kinship.mother,person) - end - if kinship.partner != nothing - resetPartner!(kinship.partner) - partner.partner = person - end - if length(kinship.children) > 0 - for child in kinship.children - setAsParentChild!(person,child) - end - end - person - end # Person Cor +@composite @kwdef mutable struct Person + + BasicInfo... + Kinship{Person}... + Maternity... + Work... + Care... + Class... + Dependency{Person}... + + pos::House{Person, Town} = undefinedHouse + # undefined Person + function Person(::Nothing) + p = new() + p.age = -1 + p + end + + # default constructor + Person(args...) = new(args...) end # struct Person + # delegate functions to components # and export accessors @delegate_onefield Person pos [getHomeTown, getHomeTownName] -@export_forward Person info [age, gender, alive] -@delegate_onefield Person info [isFemale, isMale, agestep!, agestepAlive!, hasBirthday, yearsold] - -@export_forward Person kinship [father, mother, partner, children] -@delegate_onefield Person kinship [hasChildren, addChild!, isSingle, parents, siblings] - -@delegate_onefield Person maternity [startMaternity!, stepMaternity!, endMaternity!, - isInMaternity, maternityDuration] - -@export_forward Person work [status, outOfTownStudent, newEntrant, initialIncome, finalIncome, +@export_forward Person [age, gender, alive] +@export_forward Person [father, mother, partner, children, pTime] +@export_forward Person [status, outOfTownStudent, newEntrant, initialIncome, finalIncome, wage, income, potentialIncome, jobTenure, schedule, workingHours, weeklyTime, - availableWorkingHours, workingPeriods, pension] -@delegate_onefield Person work [setEmptyJobSchedule!, setFullWeeklyTime!] - -@export_forward Person care [careNeedLevel, socialWork, childWork] - -@export_forward Person class [classRank] -@delegate_onefield Person class [addClassRank!] - -@export_forward Person dependencies [guardians, dependents, provider, providees] -@delegate_onefield Person dependencies [isDependent, hasDependents, hasProvidees] - -"costum @show method for Agent person" -function Base.show(io::IO, person::Person) - print(person.info) - undefined(person.pos) ? nothing : print(" @ House id : $(person.pos.id)") - print(person.kinship) - println() -end - -#Base.show(io::IO, ::MIME"text/plain", person::Person) = Base.show(io,person) - -"Constructor with default values" - -Person(pos,age; gender=unknown, - father=nothing,mother=nothing, - partner=nothing,children=Person[]) = - Person(pos,BasicInfoBlock(;age, gender), - KinshipBlock{Person}(father,mother,partner,children), - MaternityBlock(false, 0), - WorkBlock(), - CareBlock(0, 0, 0), - ClassBlock(0), DependencyBlock{Person}()) - - -"Constructor with default values" -Person(;pos=undefinedHouse,age=0, - gender=unknown, - father=nothing,mother=nothing, - partner=nothing,children=Person[]) = - Person(pos,BasicInfoBlock(;age,gender), - KinshipBlock{Person}(father,mother,partner,children), - MaternityBlock(false, 0), - WorkBlock(), - CareBlock(0, 0, 0), - ClassBlock(0), DependencyBlock{Person}()) - + availableWorkingHours, workingPeriods, workExperience, pension] +@export_forward Person [careNeedLevel, socialWork, childWork] +@export_forward Person [classRank, parentClassRank] +@export_forward Person [guardians, dependents, provider, providees] const PersonHouse = House{Person, Town} +const PersonTown = Town{PersonHouse} +const undefinedTown = PersonTown((-1,-1), 0.0) const undefinedHouse = PersonHouse(undefinedTown, (-1, -1)) +const undefinedPerson = Person(nothing) + +undefined(::Type{PersonHouse}) = undefinedHouse +undefined(::Type{PersonTown}) = undefinedTown +undefined(::Type{Person}) = undefinedPerson +isUndefined(h::PersonHouse) = h == undefinedHouse +isUndefined(t::PersonTown) = t == undefinedTown +isUndefined(p::Person) = p == undefinedPerson + "associate a house to a person, removes person from previous house" function moveToHouse!(person::Person,house) - if ! undefined(person.pos) + if ! isUndefined(person.pos) removeOccupant!(person.pos, person) end @@ -164,7 +98,7 @@ end "reset house of a person (e.g. became dead)" function resetHouse!(person::Person) - if ! undefined(person.pos) + if ! isUndefined(person.pos) removeOccupant!(person.pos, person) end @@ -174,10 +108,21 @@ end livingTogether(person1, person2) = person1.pos == person2.pos +"Whether the person shares their house with a non-dependent, non-guardian. Note that this includes spouses and spouses' children." +function livesInSharedHouse(person) + for p in person.pos.occupants + if p != person && ! (p in guardians(person)) && ! (p in dependents(person)) + return true + end + end + + false +end + areParentChild(person1, person2) = person1 in children(person2) || person2 in children(person1) -areSiblings(person1, person2) = father(person1) == father(person2) != nothing || - mother(person1) == mother(person2) != nothing +areSiblings(person1, person2) = father(person1) == father(person2) != undefinedPerson || + mother(person1) == mother(person2) != undefinedPerson related1stDegree(person1, person2) = areParentChild(person1, person2) || areSiblings(person1, person2) @@ -191,8 +136,8 @@ householdIncomePerCapita(person) = householdIncome(person) / length(person.pos.o function setAsParentChild!(child::Person,parent::Person) @assert isMale(parent) || isFemale(parent) @assert age(child) < age(parent) - @assert (isMale(parent) && father(child) == nothing) || - (isFemale(parent) && mother(child) == nothing) + @assert (isMale(parent) && isUndefined(father(child))) || + (isFemale(parent) && isUndefined(mother(child))) addChild!(parent, child) setParent!(child, parent) @@ -204,9 +149,11 @@ end function resetPartner!(person) other = partner(person) - if other != nothing - partner!(person, nothing) - partner!(other, nothing) + if !isUndefined(other) + partner!(person, undefinedPerson) + pTime!(person, 0) + partner!(other, undefinedPerson) + pTime!(other, 0) end nothing end @@ -279,6 +226,9 @@ isOrphan(person) = !canLiveAlone(person) && !isDependent(person) function setAsGuardianDependent!(guardian, dependent) push!(dependents(guardian), dependent) push!(guardians(dependent), guardian) + + # set class rank to maximum of guardians' + parentClassRank!(dependent, maximum(classRank, guardians(dependent))) nothing end @@ -297,9 +247,10 @@ function resolveDependency!(guardian, dependent) error("inconsistent dependency!") end deleteat!(guards, idx_g) + nothing end - +"Dissolve all guardian-dependent relationships of `person`" function setAsIndependent!(person) if !isDependent(person) return @@ -317,12 +268,12 @@ end # then something is seriously wrong function checkConsistencyDependents(person) for guard in guardians(person) - @assert guard != nothing && alive(guard) + @assert !isUndefined(guard) && alive(guard) @assert person in dependents(guard) end for dep in dependents(person) - @assert dep != nothing && alive(dep) + @assert !isUndefined(dep) @assert age(dep) < 18 @assert person.pos == dep.pos @assert person in guardians(dep) @@ -331,7 +282,7 @@ end function setAsProviderProvidee!(prov, providee) - @assert provider(providee) == nothing + @assert isUndefined(provider(providee)) @assert !(providee in providees(prov)) push!(providees(prov), providee) provider!(providee, prov) @@ -339,13 +290,13 @@ function setAsProviderProvidee!(prov, providee) end function setAsSelfproviding!(person) - if provider(person) == nothing + if isUndefined(provider(person)) return end provs = providees(provider(person)) deleteat!(provs, findfirst(==(person), provs)) - provider!(person, nothing) + provider!(person, undefinedPerson) nothing end @@ -353,13 +304,54 @@ end function maxParentRank(person) f = father(person) m = mother(person) - if f == m == nothing + if f == m == undefinedPerson classRank(person) - elseif f == nothing + elseif f == undefinedPerson classRank(m) - elseif m == nothing + elseif m == undefinedPerson classRank(f) else max(classRank(m), classRank(f)) end end + +function Utilities.dump_header(io, p::Person, FS) + print(io, "id", FS, "house", FS) + Utilities.dump_header(io, p.info, FS); print(io, FS) + Utilities.dump_header(io, p.kinship, FS); print(io, FS) + Utilities.dump_header(io, p.maternity, FS); print(io, FS) + Utilities.dump_header(io, p.work, FS); print(io, FS) + Utilities.dump_header(io, p.care, FS); print(io, FS) + Utilities.dump_header(io, p.class, FS); print(io, FS) + Utilities.dump_header(io, p.dependencies, FS) +end + +# links to objects are dumped as object id +function Utilities.dump_property(io, prop::Person, FS="\t", ES=",") + print(io, objectid(prop)) +end + +function Utilities.dump_property(io, prop::House, FS="\t", ES=",") + print(io, objectid(prop)) +end + +function Utilities.dump_property(io, prop::Union{Person, Nothing}, FS="\t", ES=",") + if prop == nothing + print(io, 0) + else + Utilities.dump_property(io, prop, FS, ES) + end +end + +function Utilities.dump(io, person::Person, FS="\t", ES=",") + print(io, objectid(person), FS) + Utilities.dump_property(io, person.pos, FS, ES); print(io, FS) + Utilities.dump(io, person.info, FS, ES); print(io, FS) + Utilities.dump(io, person.kinship, FS, ES); print(io, FS) + Utilities.dump(io, person.maternity, FS, ES); print(io, FS) + Utilities.dump(io, person.work, FS, ES); print(io, FS) + Utilities.dump(io, person.care, FS, ES); print(io, FS) + Utilities.dump(io, person.class, FS, ES); print(io, FS) + Utilities.dump(io, person.dependencies, FS, ES) +end + diff --git a/src/agents/town.jl b/src/agents/town.jl index 1c8e370..923ab9f 100644 --- a/src/agents/town.jl +++ b/src/agents/town.jl @@ -1,6 +1,7 @@ -export Town, undefinedTown, TownLocation +export Town, TownLocation export isAdjacent8, adjacent8Towns, manhattanDistance + """ Specification of a Town agent type. @@ -15,37 +16,22 @@ Type Town to extend from AbstractAXgent. const TownLocation = NTuple{2,Int} -struct Town <: AbstractXAgent - id::Int +struct Town{H} pos::TownLocation - name::String # does not look necessary - # lha::Array{Float64,1} # local house allowance - # a vector of size 4 each corresponding the number of bed rooms density::Float64 # relative population density w.r.t. the town with the highest density - - "" - function Town(pos::TownLocation,name::String,density::Float64) - #global IDCOUNTER = IDCOUNTER + 1 - # idcounter = getIDCOUNTER() - # new(IDCOUNTER,pos,name,density) - new(getIDCOUNTER(),pos,name,density) - end - + houses :: Vector{H} + adjacent :: Vector{Town{H}} end # Town +Town{H}(pos, density) where {H} = Town(pos, density, H[], Town{H}[]) + "costum show method for Town" function Base.show(io::IO, town::Town) print(" Town $(town.id) ") - isempty(town.name) ? nothing : print(" $(town.name) ") print("@ $(town.pos)") println(" density: $(town.density)") end -# Base.show(io::IO, ::MIME"text/plain", town::Town) = Base.show(io,town) - -Town(pos;name="",density=0.0) = Town(pos,name,density) - -const undefinedTown = Town((-1,-1),"",0.0) isAdjacent8(town1, town2) = abs(town1.pos[1] - town2.pos[1]) <= 1 && diff --git a/src/agents/world.jl b/src/agents/world.jl index 63f5aea..5b47bcf 100644 --- a/src/agents/world.jl +++ b/src/agents/world.jl @@ -2,18 +2,13 @@ using Memoization export adjacent8Towns, findHousesInTown, emptyHouses, emptyHousesInTown -# memoization not really necessary for low number of towns, but why not -"Find all towns adjacent to `town` (von Neumann neighbourhood). Memoized for efficiency - empty cache when topology changes." -@memoize adjacent8Towns(town, towns) = [ t for t in towns if isAdjacent8(town, t) ] - - -"Find all houses belonging to a specific town. Memoization should make this really fast after a few iterations. Empty cache only when houses appear/disappear." -@memoize findHousesInTown(t, allHouses) = - [ house for house in allHouses if town(house) == t ] +"Find all towns adjacent to `town` (von Neumann neighbourhood)." +adjacent8Towns(town) = town.adjacent +"Find all houses belonging to a specific town." +findHousesInTown(t) = t.houses emptyHouses(allHouses) = [ house for house in allHouses if isEmpty(house) ] -emptyHousesInTown(town, allHouses) = - [ house for house in findHousesInTown(town, allHouses) if isEmpty(house) ] - +emptyHousesInTown(town) = + [ house for house in findHousesInTown(town) if isEmpty(house) ] diff --git a/src/demography/data.jl b/src/demography/data.jl new file mode 100644 index 0000000..f4b4672 --- /dev/null +++ b/src/demography/data.jl @@ -0,0 +1,29 @@ +using CSV +using Tables + + +export loadDemographyData, DemographyData + +struct DemographyData + initialAgePyramid :: Vector{Vector{Float64}} + pre51Fertility :: Matrix{Float64} + fertility :: Matrix{Float64} + pre51Deaths :: Matrix{Float64} + deathFemale :: Matrix{Float64} + deathMale :: Matrix{Float64} +end + +function loadDemographyData(apFName, pre51FertFName, fertFName, + pre51DeathsFName, deathFFName, deathMFName) + + agePyramid = CSV.File(apFName, header=0) |> Tables.matrix + pre51Fert = CSV.File(pre51FertFName, header=0) |> Tables.matrix + fert = CSV.File(fertFName, header=0) |> Tables.matrix + pre51Deaths = CSV.File(pre51DeathsFName, header=0) |> Tables.matrix + deathFemale = CSV.File(deathFFName, header=0) |> Tables.matrix + deathMale = CSV.File(deathMFName, header=0) |> Tables.matrix + + DemographyData([cumsum(agePyramid[:, 1]), cumsum(agePyramid[:, 2])], + pre51Fert, fert, pre51Deaths, deathFemale, deathMale) +end + diff --git a/src/lpm/demography/demographypars.jl b/src/demography/parameters.jl similarity index 85% rename from src/lpm/demography/demographypars.jl rename to src/demography/parameters.jl index 5c602be..f69b2aa 100644 --- a/src/lpm/demography/demographypars.jl +++ b/src/demography/parameters.jl @@ -1,5 +1,5 @@ using Parameters -export MapPars, PopulationPars, DemographyPars, DivorcePars, WorkPars +export MapPars, PopulationPars, DivorcePars, WorkPars, ModelPars "Parameters describing map properties" @with_kw mutable struct MapPars @@ -100,15 +100,10 @@ end # UKMapPars "Parameters related to population setup and dynamics" @with_kw mutable struct PopulationPars - baseDieProb::Float64 = 0.0001 - babyDieProb::Float64 = 0.005 - femaleAgeDieProb::Float64 = 0.00019 - femaleAgeScaling::Float64 = 15.5 - femaleMortalityBias::Float64 = 0.85 + startTime :: Rational{Int} = 1920 + finishTime :: Rational{Int} = 2040 initialPop::Int = 5000 # Number of females or males in the initial population - maleAgeDieProb::Float64 = 0.00021 - maleAgeScaling::Float64 = 14.0 - maleMortalityBias::Float64 = 0.8 + initialPMales :: Float64 = 0.477 # from 1921 census # a population of males to be randomly generated in the # range of minStartAge - maxStartAge maxStartAge::Int = 45 @@ -118,16 +113,27 @@ end # UKMapPars startProbMarried::Float64 = 0.8 startProbOrphan::Float64 = 0.01 - cumProbClasses::Vector{Float64} = cumsum([0.2, 0.35, 0.25, 0.15, 0.05]) + baseDieProb::Float64 = 0.0001 + babyDieProb::Float64 = 0.005 + femaleAgeDieProb::Float64 = 0.00019 + femaleAgeScaling::Float64 = 15.5 + femaleMortalityBias::Float64 = 0.85 + maleAgeDieProb::Float64 = 0.00021 + maleAgeScaling::Float64 = 14.0 + maleMortalityBias::Float64 = 0.8 + + cumProbClasses::Vector{Float64} = cumsum([0.2, 0.2, 0.2, 0.2, 0.2]) +# cumProbClasses::Vector{Float64} = cumsum([0.2, 0.35, 0.25, 0.15, 0.05]) end # PopulationPars "Parameters related to reproduction" @with_kw mutable struct BirthPars fertilityBias::Float64 = 0.9 + prevChildFertBias::Float64 = 0.9 growingPopBirthProb::Float64 = 0.215 - maxPregnancyAge::Int = 42 - minPregnancyAge::Int = 17 + maxPregnancyAge::Int = 50 + minPregnancyAge::Int = 16 end @@ -144,6 +150,7 @@ end incomeFinalLevels :: Vector{Float64} = [12.0, 16.0, 25.0, 40.0, 60.0] incomeGrowthRate :: Vector{Float64} = [0.4/12.0, 0.35/12.0, 0.3/12.0, 0.25/12.0, 0.2/12.0] workingAge :: Vector{Int} = [16, 18, 20, 22, 24] + "working hours by care requirement" weeklyHours :: Vector{Float64} = [40.0, 20.0, 10.0, 0.0, 0.0] constantIncomeParam :: Float64 = 50.0 constantEduParam :: Float64 = 4.0 @@ -151,6 +158,7 @@ end eduRankSensitivity :: Float64 = 4.0 careEducationParam :: Float64 = 0.0 workDiscountingTime :: Float64 = 1.0 + moveOutProb :: Float64 = 0.1 end @@ -180,30 +188,57 @@ end betaSocExp :: Float64 = 2.0 rankGenderBias :: Float64 = 0.5 "prob dist of age difference" - deltaAgeProb :: Vector{Float64} = [0.0, 0.1, 0.25, 0.4, 0.2, 0.05] + modeAgeDiff :: Float64 = 2.0 + maleOlderFactor :: Float64 = 0.1 + maleYoungerFactor :: Float64 = 0.3 bridesChildrenExp :: Float64 = 0.5 end +"Social care" +@with_kw mutable struct CarePars + careBias :: Float64 = 0.9 + femaleAgeCareScaling :: Float64 = 19.0 + maleAgeCareScaling :: Float64 = 18.0 + personCareProb :: Float64 = 0.0008 + baseCareProb :: Float64 = 0.0002 + "care demand dependent on need level" + careDemandInHours :: Vector{Int} = [ 0, 14, 28, 56, 84 ] + careTransitionRate :: Float64 = 0.7 + zeroYearCare :: Int = 80 + childCareDemand :: Int = 168 + freeChildCareHoursPreSchool :: Int = 24 + freeChildCareHoursSchool :: Int = 32 + "weekly care supply for child, teenager, student, worker, retired, unemployed" + careSupplyByStatus :: Vector{Int} = [ 0, 10, 24, 32, 60, 48 ] + careQuantum :: Int = 2 +end + "Data files" @with_kw mutable struct DataPars datadir :: String = "data" + iniAgeFName :: String = "age_pyramid_1921_5y.csv" + pre51FertFName :: String = "birthrate_early.csv" fertFName :: String = "babyrate.txt.csv" + pre51DeathsFName :: String = "deathrates_early.csv" deathFFName :: String = "deathrate.fem.csv" deathMFName :: String = "deathrate.male.csv" end -struct DemographyPars + + +struct ModelPars mappars :: MapPars poppars :: PopulationPars birthpars :: BirthPars workpars :: WorkPars divorcepars :: DivorcePars marriagepars :: MarriagePars + carepars :: CarePars datapars :: DataPars end -DemographyPars() = DemographyPars(MapPars(), PopulationPars(), BirthPars(), WorkPars(), - DivorcePars(), MarriagePars(), DataPars()) +ModelPars() = ModelPars(MapPars(), PopulationPars(), BirthPars(), WorkPars(), + DivorcePars(), MarriagePars(), CarePars(), DataPars()) diff --git a/src/demography/setup/map.jl b/src/demography/setup/map.jl new file mode 100644 index 0000000..141d560 --- /dev/null +++ b/src/demography/setup/map.jl @@ -0,0 +1,53 @@ + +export createTowns, initializeHousesInTowns + + +function createTowns(pars) + towns = Matrix{PersonTown}(undef, pars.mapGridXDimension, pars.mapGridYDimension) + + for y in 1:pars.mapGridYDimension, x in 1:pars.mapGridXDimension + towns[x, y] = PersonTown((x,y), pars.map[y,x]) + end + + for t in towns + x = t.pos[1] + y = t.pos[2] + for xx in x-1:x+1, yy in y-1:y+1 + if xx == x && yy == y + continue + end + + if xx < 1 || xx > size(towns, 1) || yy < 1 || yy > size(towns, 2) + continue + end + + push!(t.adjacent, towns[xx, yy]) + end + +# previous version produced the same neighbours, but in a different order +# adj2 = [t2 for t2 in towns if isAdjacent8(t, t2) ] +# @assert issetequal(adj2, t.adjacent) + end + + vec(towns) +end + + +"initialize houses in a given set of towns" +function initializeHousesInTowns(towns::Vector{PersonTown}, pars) + houses = PersonHouse[] + + for town in towns + adjustedDensity = town.density * pars.mapDensityModifier + + for hx in 1:pars.townGridDimension, hy in 1:pars.townGridDimension + if rand() < adjustedDensity + house = PersonHouse(town, (hx,hy)) + push!(houses,house) + push!(town.houses, house) + end + end # for hx + end # for town + + houses +end # function initializeHousesInTwons diff --git a/src/demography/setup/mapPop.jl b/src/demography/setup/mapPop.jl new file mode 100644 index 0000000..b59de52 --- /dev/null +++ b/src/demography/setup/mapPop.jl @@ -0,0 +1,31 @@ +using Random + +export assignCouplesToHouses! + +"Randomly assign a population of couples to non-inhebted set of houses" +function assignCouplesToHouses!(population::Array{Person}, houses::Array{PersonHouse}) + women = [ person for person in population if isFemale(person) ] + + randomhouses = shuffle(houses) + + for woman in women + house = pop!(randomhouses) + + moveToHouse!(woman, house) + if !isSingle(woman) + moveToHouse!(partner(woman), house) + end + + for child in dependents(woman) + moveToHouse!(child, house) + end + end # for person + + for person in population + if person.pos == undefinedHouse + @assert isMale(person) + @assert length(randomhouses) >= 1 + moveToHouse!(person, pop!(randomhouses)) + end + end +end # function assignCouplesToHouses diff --git a/src/lpm/demography/Create.jl b/src/demography/setup/population.jl similarity index 62% rename from src/lpm/demography/Create.jl rename to src/demography/setup/population.jl index 278d8ab..fdc1c20 100644 --- a/src/lpm/demography/Create.jl +++ b/src/demography/setup/population.jl @@ -1,30 +1,3 @@ -""" - -""" - -module Create - -using Distributions - -using Utilities -using XAgents - -export createTowns, createPopulation, createPyramidPopulation -### - -function createTowns(pars) - - uktowns = Town[] - - for y in 1:pars.mapGridYDimension - for x in 1:pars.mapGridXDimension - town = Town((x,y),density=pars.map[y,x]) - push!(uktowns,town) - end - end - - uktowns -end # return agents with age in interval minAge, maxAge # assumes pop is sorted by age @@ -51,26 +24,36 @@ function ageInterval(pop, minAge, maxAge) idx_start, idx_end end +function randAge(pyramid, gender) + data = pyramid[gender == male ? 1 : 2] + r = rand() * data[end] + i = searchsortedfirst(data, r) + mi = (i-1) * 5 * 12 # we are working in months + ma = i * 5 * 12 - 1 + rand(mi:ma) // 12 +end + -function createPyramidPopulation(pars) +function createPyramidPopulation(pars, pyramid) population = Person[] men = Person[] women = Person[] # age pyramid - dist = TriangularDist(0, pars.maxStartAge * 12, 0) + #dist = TriangularDist(0, pars.maxStartAge * 12, 0) for i in 1:pars.initialPop # surplus of babies and toddlers, lower bit of age pyramid - if i < pars.startBabySurplus - age = rand(1:36) // 12 - else - age = floor(Int, rand(dist)) // 12 - end + #if i < pars.startBabySurplus + # age = rand(1:36) // 12 + #else + # age = floor(Int, rand(dist)) // 12 + #end - gender = Bool(rand(0:1)) ? male : female + gender = rand() < pars.initialPMales ? male : female + age = randAge(pyramid, gender) - person = Person(undefinedHouse, age; gender) + person = Person(;age, gender) if age < 18 push!(population, person) else @@ -149,8 +132,8 @@ function createPyramidPopulation(pars) population end -function createPopulation(pars) +function createUniformPopulation(pars) population = Person[] for i in 1 : pars.initialPop @@ -168,8 +151,8 @@ function createPopulation(pars) # birthYear = properties[:startYear] - ageMale/Female # birthMonth = rand((1:12)) - newMan = Person(undefinedHouse,rageMale,gender=male) - newWoman = Person(undefinedHouse,rageFemale,gender=female) + newMan = Person(age=rageMale, gender=male) + newWoman = Person(age=rageFemale, gender=female) setAsPartners!(newMan,newWoman) push!(population,newMan); push!(population,newWoman) @@ -178,7 +161,62 @@ function createPopulation(pars) population -end # createPopulation +end # createUniformPopulation + + +function initClass!(person, pars) + p = rand() + class = searchsortedfirst(pars.cumProbClasses, p)-1 + classRank!(person, class) + nothing +end + + +function initWork!(person, pars) + if age(person) < pars.ageTeenagers + status!(person, WorkStatus.child) + return + end + if age(person) < pars.ageOfAdulthood + status!(person, WorkStatus.teenager) + return + end + if age(person) >= pars.ageOfRetirement + status!(person, WorkStatus.retired) + return + end + + class = classRank(person)+1 + + if age(person) < pars.workingAge[class] + status!(person, WorkStatus.student) + return + end -end # module Create + status!(person, WorkStatus.worker) + + workingTime = 0 + for i in pars.workingAge[class]:floor(Int, age(person)) + workingTime *= pars.workDiscountingTime + workingTime += 1 + end + + workExperience!(person, workingTime) + workingPeriods!(person, workingTime) + + dKi = rand(Normal(0, pars.wageVar)) + initialWage = pars.incomeInitialLevels[class] * exp(dKi) + dKf = rand(Normal(dKi, pars.wageVar)) + finalWage = pars.incomeFinalLevels[class] * exp(dKf) + + initialIncome!(person, initialWage) + finalIncome!(person, finalWage) + + wage!(person, computeWage(person, pars)) + income!(person, wage(person) * pars.weeklyHours[careNeedLevel(person)+1]) + potentialIncome!(person, income(person)) + jobTenure!(person, rand(1:50)) + + nothing +end diff --git a/src/lpm/demography/simulate/ageTransition.jl b/src/demography/simulate/ageTransition.jl similarity index 60% rename from src/lpm/demography/simulate/ageTransition.jl rename to src/demography/simulate/ageTransition.jl index bc0e6fb..04b63a0 100644 --- a/src/lpm/demography/simulate/ageTransition.jl +++ b/src/demography/simulate/ageTransition.jl @@ -1,7 +1,5 @@ using Distributions -using XAgents - export selectAgeTransition, ageTransition!, selectWorkTransition, workTransition! @@ -24,6 +22,11 @@ function ageTransition!(person, time, model, pars) # person.yearInTown += 1 #end agestep!(person) + + # TODO parameterise dt + if !isSingle(person) + pTime!(person, pTime(person) + 1//12) + end if age(person) == 18 # also updates guardian @@ -34,6 +37,21 @@ end selectWorkTransition(person, pars) = alive(person) && status(person) != WorkStatus.retired && hasBirthday(person) +function computeWage(person, pars) + # original formula + # c = log(I/F) + # wage = F * exp(c * exp(-1 * r * e)) + + fI = finalIncome(person) + iI = initialIncome(person) + + wage = fI * (iI/fI)^exp(-1 * pars.incomeGrowthRate[classRank(person)+1] * workExperience(person)) + + dK = rand(Normal(0, pars.wageVar)) + + wage * exp(dK) +end + function workTransition!(person, time, model, pars) if age(person) == pars.ageTeenagers @@ -43,7 +61,7 @@ function workTransition!(person, time, model, pars) if age(person) == pars.ageOfAdulthood status!(person, WorkStatus.student) - #class!(person, 0) + classRank!(person, 0) if rand() < pars.probOutOfTownStudent outOfTownStudent!(person, true) @@ -61,5 +79,17 @@ function workTransition!(person, time, model, pars) dK = rand(Normal(0, pars.wageVar)) pension!(person, shareWorkingTime * exp(dK)) + return + end + + if status(person) == WorkStatus.worker && !isInMaternity(person) + # we assume full work load at this point + # in original: availableWorkingHours/workingHours + workingPeriods!(person, workingPeriods(person)+1) + # in original: availableWorkingHours/pars.weeklyHours[0] + workExperience!(person, workExperience(person)+1) + wage!(person, computeWage(person, pars)) + # no care, therefore full time + income!(person, wage(person) * pars.weeklyHours[careNeedLevel(person)+1]) end end diff --git a/src/lpm/demography/simulate/allocate.jl b/src/demography/simulate/allocate.jl similarity index 57% rename from src/lpm/demography/simulate/allocate.jl rename to src/demography/simulate/allocate.jl index 752ec0a..4278596 100644 --- a/src/lpm/demography/simulate/allocate.jl +++ b/src/demography/simulate/allocate.jl @@ -1,32 +1,21 @@ - -#= -An initial design for findNewHouse*(*) interfaces (subject to incremental - modification, simplifcation and tuning) -=# - -using Memoization - -using XAgents - export findEmptyHouseInTown, findEmptyHouseInOrdAdjacentTown, findEmptyHouseAnywhere, movePeopleToEmptyHouse!, movePeopleToHouse! function selectHouse(list) if isempty(list) - return nothing + return undefinedHouse end rand(list) end -findEmptyHouseInTown(town, allHouses) = selectHouse(emptyHousesInTown(town, allHouses)) - +findEmptyHouseInTown(town) = selectHouse(emptyHousesInTown(town)) function findEmptyHouseInOrdAdjacentTown(town, allHouses, allTowns) - adjTowns = adjacent8Towns(town, allTowns) - emptyHouses = [ house for town in adjTowns - for house in findHousesInTown(town, allHouses) if isEmpty(house) ] + adjTowns = [town; adjacent8Towns(town)] + emptyHouses = [ house for t in adjTowns + for house in findHousesInTown(t) if isEmpty(house) ] selectHouse(emptyHouses) end @@ -47,22 +36,22 @@ function movePeopleToHouse!(people, house) end # people[1] determines centre of search radius -function movePeopleToEmptyHouse!(people, dmax, allHouses, allTowns=Town[]) - newhouse = nothing +function movePeopleToEmptyHouse!(people, dmax, allHouses, allTowns=PersonTown[]) + newhouse = undefinedHouse if dmax == :here - newhouse = findEmptyHouseInTown(people[1].pos,allHouses) + newhouse = findEmptyHouseInTown(getHomeTown(people[1])) end - if dmax == :near || newhouse == nothing - newhouse = findEmptyHouseInOrdAdjacentTown(people[1].pos,allHouses,allTowns) + if dmax == :near || newhouse == undefinedHouse + newhouse = findEmptyHouseInOrdAdjacentTown(getHomeTown(people[1]), allHouses, allTowns) end - if dmax == :far || newhouse == nothing + if dmax == :far || newhouse == undefinedHouse newhouse = findEmptyHouseAnywhere(allHouses) end - if newhouse != nothing + if newhouse != undefinedHouse movePeopleToHouse!(people, newhouse) return newhouse end - nothing + undefinedHouse end diff --git a/src/demography/simulate/birth.jl b/src/demography/simulate/birth.jl new file mode 100644 index 0000000..ca77ae4 --- /dev/null +++ b/src/demography/simulate/birth.jl @@ -0,0 +1,189 @@ +using Utilities + +export selectBirth, birth! + +isFertileWoman(p, pars) = isFemale(p) && pars.minPregnancyAge <= age(p) <= pars.maxPregnancyAge +canBePregnant(p) = !isSingle(p) && ageYoungestAliveChild(p) > 1 +isPotentialMother(p, pars) = isFertileWoman(p, pars) && canBePregnant(p) + +mutable struct BirthCache{PERSON} + potentialMothers :: Vector{PERSON} + pPotentialMotherInFertWAndAge :: Vector{Float64} + classBias :: Vector{Float64} + nChBias :: Matrix{Float64} +end + +BirthCache{T}() where {T} = BirthCache(T[], Float64[], Float64[], zeros(5, 5)) + +function birthPreCalc!(model, pars) + pc = model.birthCache + empty!(pc.potentialMothers) + for a in model.pop + if isPotentialMother(a, pars) + push!(pc.potentialMothers, a) + end + end + + resize!(pc.pPotentialMotherInFertWAndAge, 150) + fill!(pc.pPotentialMotherInFertWAndAge, 0) + cbp = zeros(Int, 150) + for a in model.pop + if isFertileWoman(a, pars) + years = yearsold(a) + cbp[years] += 1 + if canBePregnant(a) + pc.pPotentialMotherInFertWAndAge[years] += 1 + end + end + end + for (i, n) in enumerate(cbp) + pc.pPotentialMotherInFertWAndAge[i] /= n > 0 ? n : Inf + end + + nPerClass = zeros(5) + for p in pc.potentialMothers + nPerClass[classRank(p)+1] += 1 + end + + pcpm = copy(nPerClass) + if length(pc.potentialMothers) > 0 + pcpm ./= length(pc.potentialMothers) + end + resize!(pc.classBias, 5) + preCalcRateBias!(c->pcpm[c+1], 0:4, pars.fertilityBias, pc.classBias, 1) + + pncpmc = zeros(5, 5) + for p in pc.potentialMothers + c = classRank(p) + nc = min(4, nChildren(p)) + pncpmc[c+1, nc+1] += 1 + end + + fill!(pc.nChBias, 0.0) + for class in 0:4 + for n in 0:4 + pncpmc[class+1, n+1] /= nPerClass[class+1] + end + preCalcRateBias!(n -> pncpmc[class+1, n+1], 0:4, pars.prevChildFertBias, + view(pc.nChBias, class+1, :), 1) + end +end + +"Proportion of women that can get pregnant in entire population." +function pPotentialMotherInAllPop(model, pars) + n = length(model.birthCache.potentialMothers) + + n / length(model.pop) +end + + +function computeBirthProb(woman, parameters, model, currstep) + (curryear,currmonth) = date2yearsmonths(currstep) + currmonth = currmonth + 1 # adjusting 0:11 => 1:12 + + womanRank = classRank(woman) + if status(woman) == WorkStatus.student + womanRank = parentClassRank(woman) + end + + ageYears = yearsold(woman) + fertAge = ageYears-parameters.minPregnancyAge+1 + + if curryear < 1951 + # number of children per uk resident and year + rawRate = model.pre51Fertility[Int(curryear-parameters.startTime+1)] / + # scale by number of women that can actually get pregnant + pPotentialMotherInAllPop(model, parameters) * + # and multiply with age-specific fertility factor + model.fertFByAge51[fertAge] + else + # fertility rates are stored as P(pregnant) per year and age + rawRate = model.fertility[fertAge, curryear-1950] / + model.birthCache.pPotentialMotherInFertWAndAge[ageYears] + end + + # fertility bias by class + birthProb = rawRate * model.birthCache.classBias[womanRank+1] + + # fertility bias by number of previous children + birthProb *= model.birthCache.nChBias[womanRank+1, min(4,nChildren(woman))+1] + + min(1.0, birthProb) +end # computeBirthProb + + +function effectsOfMaternity!(woman, pars) + startMaternity!(woman) + + workingHours!(woman, 0) + income!(woman, 0) + potentialIncome!(woman, 0) + availableWorkingHours!(woman, 0) + # commented in sim.py: + # woman.weeklyTime = [[0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12] + # sets all weeklyTime slots to 1 + # TODO copied from the python code, but does it make sense? + setFullWeeklyTime!(woman) + #= TODO + woman.maxWeeklySupplies = [0, 0, 0, 0] + woman.residualDailySupplies = [0]*7 + woman.residualWeeklySupplies = [x for x in woman.maxWeeklySupplies] + =# + + # TODO not necessarily true in many cases + if provider(woman) == undefinedPerson + setAsProviderProvidee!(partner(woman), woman) + end + + nothing +end + + +selectBirth(person, parameters) = isFertileWoman(person, parameters) && !isSingle(person) && + ageYoungestAliveChild(person) > 1 + + +function birth!(woman, currstep, model, parameters, addBaby!) + birthProb = computeBirthProb(woman, parameters, model, currstep) + + assumption() do + @assert isFemale(woman) + @assert ageYoungestAliveChild(woman) > 1 + @assert !isSingle(woman) + @assert age(woman) >= parameters.minPregnancyAge + @assert age(woman) <= parameters.maxPregnancyAge + @assert birthProb >= 0 + end + + #= + The following code is commented in the python code: + #baseRate = self.baseRate(self.socialClassShares, self.p['fertilityBias'], rawRate) + #fertilityCorrector = (self.socialClassShares[woman.classRank] - self.p['initialClassShares'][woman.classRank])/self.p['initialClassShares'][woman.classRank] + #baseRate *= 1/math.exp(self.p['fertilityCorrector']*fertilityCorrector) + #birthProb = baseRate*math.pow(self.p['fertilityBias'], woman.classRank) + =# + + if rand() < p_yearly2monthly(limit(0.0, birthProb, 1.0)) + + baby = Person(gender=rand([male,female])) + moveToHouse!(baby, woman.pos) + setAsParentChild!(baby, woman) + if !isSingle(woman) # currently not an option + setAsParentChild!(baby, partner(woman)) + end + + # this goes first, so that we know material circumstances + effectsOfMaternity!(woman, parameters) + + setAsGuardianDependent!(woman, baby) + if !isSingle(woman) # currently not an option + setAsGuardianDependent!(partner(woman), baby) + end + setAsProviderProvidee!(woman, baby) + + addBaby!(model, baby) + end # if rand() + + nothing +end + diff --git a/src/demography/simulate/care.jl b/src/demography/simulate/care.jl new file mode 100644 index 0000000..6f3ba34 --- /dev/null +++ b/src/demography/simulate/care.jl @@ -0,0 +1,147 @@ +include("kinship.jl") + +export socialCareSupply, socialCareDemand, householdSocialCareNeed + +socialCareDemand(person, pars) = pars.careDemandInHours[careNeedLevel(person)+1] + +numCareLevels(pars) = length(pars.careDemandInHours) + + +function childCareNeed(child, model, pars) + ageC = age(child) + + childCare = + if ageC < 1 + pars.zeroYearCare + elseif ageC < 13 + pars.childCareDemand + else + 0 + end + + socCare = socialCareDemand(child, pars) + + childCare = max(0, childCare - socCare) + + schoolHours = + if ageC < 3 + 0 + elseif ageC <= 5 + pars.freeChildCareHoursPreSchool + else + pars.freeChildCareHoursSchool + end + + max(0, childCare - schoolHours) +end + + +function householdChildCareNeed(house, model, pars) + maxChildCare = 0 + for person in house.occupants + if age(person) >= 13 + continue + end + care = childCareNeed(person, model, pars) + @assert care >= 0 + maxChildCare = max(care, maxChildCare) + end + maxChildCare +end + + +function householdSocialCareNeed(house, model, pars) + careNeed = 0 + for person in house.occupants + careNeed += socialCareDemand(person, pars) + end + + careNeed +end + + +function socialCareSupply(person, pars) + if careNeedLevel(person) > 0 + return 0 + end + + s = Int(status(person)) + + pars.careSupplyByStatus[s+1] +end + +function householdSocialCareSupply(house, model, pars) + supply = 0 + for person in house.occupants + supply += socialCareSupply(person, pars) + end + + supply +end + + +function resetHouseholdCare!(house, model, pars) + cn = householdChildCareNeed(house, model, pars) + sn = householdSocialCareNeed(house, model, pars) + s = householdSocialCareSupply(house, model, pars) + + resetCare!(house, s - cn - sn) +end + +"Find all households that have positive net care supply" +function supplyHouseholds(model) + [house for house in model.houses if netCareSupply(house) > 0 ] +end + + +function buildSupplyDemandNetwork(model, pars) + links = Link{eltype(model.houses)}[] + suppliers = supplyHouseholds(model) + + for house in suppliers + # we remove replicate links for now + l = kinshipNetwork(house, model, pars) do h + # only look at households that require care + netCareSupply(h) < 0 + end + append!(links, l) + end + + links +end + + +function resolveCareSupply!(network, model, pars) + while !isempty(network) + l = rand(1:length(network)) + link = network[l] + + # provider doesn't have enough care left + if careBalance(link.t1) <= pars.careQuantum || + # providee is already fine + careBalance(link.t2) >= 0 + remove_unsorted!(network, l) + continue + end + + receiveCare!(link.t2, pars.careQuantum, link.t1) + provideCare!(link.t1, pars.careQuantum, link.t2) + + # provider doesn't have enough care left + if careBalance(link.t1) <= pars.careQuantum || + # providee is already fine + careBalance(link.t2) >= 0 + remove_unsorted!(network, l) + end + end +end + + +function socialCare!(model, pars) + for h in model.houses + resetHouseholdCare!(h, model, pars) + end + + network = buildSupplyDemandNetwork(model, pars) + resolveCareSupply!(network, model, pars) +end diff --git a/src/demography/simulate/death.jl b/src/demography/simulate/death.jl new file mode 100644 index 0000000..462490a --- /dev/null +++ b/src/demography/simulate/death.jl @@ -0,0 +1,172 @@ +using Utilities + +export death!, setDead! + +function deathProbability(baseRate, person, model, parameters) + cRank = classRank(person) + if status(person) == WorkStatus.child || status(person) == WorkStatus.student + cRank = parentClassRank(person) + end + + if isMale(person) + mortalityBias = parameters.maleMortalityBias + else + mortalityBias = parameters.femaleMortalityBias + end + + a = isMale(person) ? model.deathCache.classBias_m : model.deathCache.classBias_f + #sumClassBias(c -> model.socialCache.socialClassShares[c+1], + # 0:(length(parameters.cumProbClasses)-1), + # mortalityBias) + + if a > 0 + lowClassRate = baseRate / a + classRate = lowClassRate * mortalityBias^cRank + deathProb = classRate + else + deathProb = baseRate + end + + #= + b = 0 + for i in range(int(self.p['numCareLevels'])): + b += self.careNeedShares[classRank][i]*math.pow(self.p['careNeedBias'], (self.p['numCareLevels']-1) - i) + + if b > 0: + higherNeedRate = classRate/b + deathProb = higherNeedRate*math.pow(self.p['careNeedBias'], (self.p['numCareLevels']-1) - person.careNeedLevel) # deathProb +=# + ##### Temporarily by-passing the effect of Unmet Care Need ############# + + # The following code is already commented in the python code + # a = 0 + # for x in classPop: + # a += math.pow(self.p['unmetCareNeedBias'], 1-x.averageShareUnmetNeed) + # higherUnmetNeed = (classRate*len(classPop))/a + # deathProb = higherUnmetNeed*math.pow(self.p['unmetCareNeedBias'], 1-shareUnmetNeed) + deathProb +end # function deathProb + + +function setDead!(person) + alive!(person, false) + resetHouse!(person) + if !isSingle(person) + resolvePartnership!(partner(person),person) + end + + # dead persons are no longer dependents + setAsIndependent!(person) + + # dead persons no longer have to be provided for + setAsSelfproviding!(person) + + for p in providees(person) + provider!(p, undefinedPerson) + # TODO update provision/work status + end + empty!(providees(person)) + + # dependents are being taken care of by assignGuardian! + nothing +end + +mutable struct DeathCache + avgDieProb_m :: Float64 + avgDieProb_f :: Float64 + classBias_m :: Float64 + classBias_f :: Float64 +end + +DeathCache() = DeathCache(0.0, 0.0, 0.0, 0.0) + +function deathPreCalc!(model, pars) + pc = model.deathCache + + s_m = 0.0 + s_f = 0.0 + n_m = 0 + n_f = 0 + for p in model.pop + if isMale(p) + s_m += ageDieProb(pars, yearsold(p), true) + n_m += 1 + else + s_f += ageDieProb(pars, yearsold(p), false) + n_f += 1 + end + end + + pc.avgDieProb_m = s_m / n_m + pc.avgDieProb_f = s_f / n_f + + pc.classBias_m = sumClassBias(c -> model.socialCache.socialClassShares[c+1], + 0:(length(pars.cumProbClasses)-1), + pars.maleMortalityBias) + pc.classBias_f = sumClassBias(c -> model.socialCache.socialClassShares[c+1], + 0:(length(pars.cumProbClasses)-1), + pars.femaleMortalityBias) +end + + +ageDieProb(pars, agep, malep) = pars.baseDieProb + (malep ? + exp(agep / pars.maleAgeScaling) * pars.maleAgeDieProb : + exp(agep / pars.femaleAgeScaling) * pars.femaleAgeDieProb) + + +# currently leaves dead agents in population +function death!(person, currstep, model, parameters) + + (curryear,currmonth) = date2yearsmonths(currstep) + currmonth += 1 # adjusting 0:11 => 1:12 + + agep = trunc(Int, age(person)) + + assumption() do + @assert alive(person) + @assert isMale(person) || isFemale(person) # Assumption + end + + if curryear < 1950 # made-up probabilities + yearIdx = trunc(Int, curryear - parameters.startTime + 1) + + if agep < 1 + rawRate = model.pre51Deaths[yearIdx, 2] / 1000.0 # infant mortality is per 1k + else + rawRate = model.pre51Deaths[yearIdx, 1] * + ageDieProb(parameters, agep, isMale(person)) / + (isMale(person) ? + model.deathCache.avgDieProb_m : model.deathCache.avgDieProb_f) + end + else + + agep = min(agep, 109) + rawRate = isMale(person) ? model.deathMale[agep+1,curryear-1950+1] : + model.deathFemale[agep+1,curryear-1950+1] + + end # currYear < 1950 + + #= + Not realized yet + classPop = [x for x in self.pop.livingPeople + if x.careNeedLevel == person.careNeedLevel] + Classes to be considered in a different module + =# + + deathProb = limit(0.0, deathProbability(rawRate, person, model, parameters), 1.0) + + #= + The following is uncommented code in the original code < 1950 + #### Temporarily by-passing the effect of unmet care need ###### + # dieProb = self.deathProb_UCN(rawRate, person.parentsClassRank, person.careNeedLevel, person.averageShareUnmetNeed, classPop) + =# + + if rand() < p_yearly2monthly(deathProb) + setDead!(person) + return true + # person.deadYear = self.year + end # rand + + false +end + diff --git a/src/lpm/demography/simulate/dependencies.jl b/src/demography/simulate/dependencies.jl similarity index 92% rename from src/lpm/demography/simulate/dependencies.jl rename to src/demography/simulate/dependencies.jl index 42e41d0..47cd443 100644 --- a/src/lpm/demography/simulate/dependencies.jl +++ b/src/demography/simulate/dependencies.jl @@ -18,7 +18,7 @@ selectAssignGuardian(person) = alive(person) && !canLiveAlone(person) && function assignGuardian!(person, time, model, pars) guard = findFamilyGuardian(person) - if guard == nothing + if isUndefined(guard) guard = findOtherGuardian(person, model.pop, pars) end @@ -27,7 +27,7 @@ function assignGuardian!(person, time, model, pars) # that are now excluded due to age won't get a chance again in the future empty!(guardians(person)) - if guard == nothing + if isUndefined(guard) return false end @@ -52,7 +52,7 @@ function findFamilyGuardian(person) # relatives of biological parents # any of these might already be guardians, but in that case they will be dead for p in pparents - if p == nothing + if isUndefined(p) continue end append!(potGuardians, parents(p)) @@ -68,13 +68,13 @@ function findFamilyGuardian(person) # potentially lots of redundancy, but just take the first # candidate that works for g in potGuardians - if g == nothing || !alive(g) || age(g) < 18 + if isUndefined(g) || !alive(g) || age(g) < 18 continue end return g end - return nothing + return undefinedPerson end function findOtherGuardian(person, people, pars) @@ -86,7 +86,7 @@ function findOtherGuardian(person, people, pars) return rand(candidates) end - return nothing + return undefinedPerson end diff --git a/src/lpm/demography/simulate/divorce.jl b/src/demography/simulate/divorce.jl similarity index 51% rename from src/lpm/demography/simulate/divorce.jl rename to src/demography/simulate/divorce.jl index 4032990..cf9f585 100644 --- a/src/lpm/demography/simulate/divorce.jl +++ b/src/demography/simulate/divorce.jl @@ -1,25 +1,31 @@ -export doDivorces!, selectDivorce, divorce! - -function divorceProbability(rawRate, pars) # ,classRank) - #= - def computeSplitProb(self, rawRate, classRank): - a = 0 - for i in range(int(self.p['numberClasses'])): - a += self.socialClassShares[i]*math.pow(self.p['divorceBias'], i) - baseRate = rawRate/a - splitProb = baseRate*math.pow(self.p['divorceBias'], classRank) - return splitProb - =# - rawRate * pars.divorceBias -end +export selectDivorce, divorce! -function divorce!(man, time, model, parameters) - applyDivorce!(man, time, model.houses, model.towns, parameters) + +mutable struct DivorceCache + classBias :: Vector{Float64} end +DivorceCache() = DivorceCache([]) -function applyDivorce!(man, time, allHouses, allTowns, parameters) - + +function divorcePreCalc!(model, pars) + pc = model.divorceCache + resize!(pc.classBias, 5) + classes = 0:(length(pars.cumProbClasses)-1) + preCalcRateBias!(classes, pars.divorceBias, pc.classBias, 1) do c + model.socialCache.socialClassShares[c+1] + end +end + + +function divorceProbability(rawRate, classRank, model, pars) + rawRate * model.divorceCache.classBias[classRank+1]#=rateBias(0:(length(pars.cumProbClasses)-1), pars.divorceBias, classRank) do c + model.socialCache.socialClassShares[c+1] + end=# +end + + +function divorce!(man, time, model, parameters) agem = age(man) assumption() do @assert isMale(man) @@ -36,9 +42,9 @@ function applyDivorce!(man, time, allHouses, allTowns, parameters) rawRate = parameters.variableDivorce * parameters.divorceModifierByDecade[ceil(Int, agem / 10)] end - divorceProb = divorceProbability(rawRate, parameters) # TODO , man.classRank) + divorceProb = divorceProbability(rawRate, classRank(man), model, parameters) - if rand() < p_yearly2monthly(divorceProb) + if rand() < p_yearly2monthly(limit(0.0, divorceProb, 1.0)) wife = partner(man) resolvePartnership!(man, wife) @@ -64,7 +70,7 @@ function applyDivorce!(man, time, allHouses, allTowns, parameters) end end # for - movePeopleToEmptyHouse!(peopleToMove, rand([:near, :far]), allHouses, allTowns) + movePeopleToEmptyHouse!(peopleToMove, rand([:near, :far]), model.houses, model.towns) return true end @@ -76,35 +82,3 @@ end selectDivorce(person, pars) = alive(person) && isMale(person) && !isSingle(person) -function doDivorces!(people, time, allHouses, allTowns, parameters) - - marriedMen = [ man for man in people if selectDivorce(man, parameters) ] - - divorced = Person[] - - for man in marriedMen - wife = partner(man) - if applyDivorce!(man, time, allHouses, allTowns, parameters) - append!(divorced,[man, wife]) - end - end - - delayedVerbose() do - println("# of divorced in current iteration $(length(divorced))") - end - - divorced -end - - -#= -Atiyah: -Draft suggestion: if an API with model argument is needed (does not seem to be the case) -an API may look like -doDivorces!(model,time,parameters) - -doDivorces!(model,time,parameters) = - doDivorces(people,time,model.houses,model.towns,parameters) - # better (to work for me as well) - # doDivorces!(population(model),time,houses(model),towns(model),divorcePars(model)) -=# diff --git a/src/demography/simulate/kinship.jl b/src/demography/simulate/kinship.jl new file mode 100644 index 0000000..93452e6 --- /dev/null +++ b/src/demography/simulate/kinship.jl @@ -0,0 +1,38 @@ +struct Link{T} + t1 :: T + t2 :: T +end + + +function kinshipNetwork(filter, house, model, pars) + LinkT = Link{typeof(house)} + conns = Vector{LinkT}() + + function checkAndAdd!(h) + if h != house && filter(h) + push!(conns, LinkT(house, h)) + end + end + + for person in house.occupants + if !isSingle(person) + checkAndAdd!(partner(person).pos) + end + + for child in children(person) + checkAndAdd!(child.pos) + end + + f = father(person) + if !isUndefined(f) + checkAndAdd!(f.pos) + end + + m = mother(person) + if !isUndefined(m) + checkAndAdd!(m.pos) + end + end + + conns +end diff --git a/src/lpm/demography/simulate/marriages.jl b/src/demography/simulate/marriages.jl similarity index 57% rename from src/lpm/demography/simulate/marriages.jl rename to src/demography/simulate/marriages.jl index e1b5309..e60e985 100644 --- a/src/lpm/demography/simulate/marriages.jl +++ b/src/demography/simulate/marriages.jl @@ -1,61 +1,62 @@ -export resetCacheMarriages, marriage!, selectMarriage +export marriage!, selectMarriage -using XAgents +using Utilities ageClass(person) = trunc(Int, age(person)/10) -@memoize Dict function shareMenNoChildren(model, ageclass :: Int) - nAll = 0 - nNoC = 0 +mutable struct MarriageCache{PERSON} + shareMenNoChildren :: Vector{Float64} + eligibleWomen :: Vector{PERSON} + weights :: Vector{Float64} +end + +MarriageCache{Person}() where {Person} = MarriageCache(Float64[], Person[], Float64[]) - for p in Iterators.filter(x->isMale(x) && ageClass(x) == ageclass, model.pop) - nAll += 1 +function marriagePreCalc!(model, pars) + pc = model.marriageCache + + resize!(pc.shareMenNoChildren, 20) + fill!(pc.shareMenNoChildren, 0.0) + nAll = zeros(Int, 20) + for p in Iterators.filter(isMale, model.pop) + ac = ageClass(p) + nAll[ac+1] += 1 # only looks at legally dependent persons (which usually are underage and # living in the same household) if !hasDependents(p) - nNoC += 1 + pc.shareMenNoChildren[ac+1] += 1 + end + end + pc.shareMenNoChildren ./= nAll + + empty!(pc.eligibleWomen) + for f in model.pop + if isFemale(f) && isSingle(f) && age(f) > pars.minPregnancyAge + push!(pc.eligibleWomen, f) end end - - nNoC / nAll -end - - -@memoize eligibleWomen(model, pars) = [f for f in model.pop if isFemale(f) && alive(f) && - isSingle(f) && age(f) > pars.minPregnancyAge] - -# reset memoization caches -# needs to be done on every time step -function resetCacheMarriages() - Memoization.empty_cache!(shareMenNoChildren) - Memoization.empty_cache!(eligibleWomen) end -function deltaAge(delta) - if delta <= -10 - 0 - elseif -10 < delta < -2 - 1 - elseif -2 < delta < 1 - 2 - elseif 1 < delta < 5 - 3 - elseif 5 < delta < 10 - 4 - else - 5 - end +function ageFactor(agem, agew, pars) + diff = Float64(agem - agew) - pars.modeAgeDiff + diff > 0 ? + 1/exp(pars.maleOlderFactor * diff^2) : + 1/exp(pars.maleYoungerFactor * diff^2) end function marryWeight(man, woman, pars) + if livingTogether(man, woman) || related1stDegree(man, woman) + return 0.0 + end + geoFactor = 1/exp(pars.betaGeoExp * geoDistance(man, woman, pars)) if status(woman) == WorkStatus.student studentFactor = pars.studentFactorParam - womanRank = maxParentRank(woman) + womanRank = parentClassRank(woman) else studentFactor = 1.0 womanRank = classRank(woman) @@ -67,71 +68,72 @@ function marryWeight(man, woman, pars) socFactor = 1/exp(betaExponent * statusDistance) - ageFactor = pars.deltaAgeProb[deltaAge(age(man) - age(woman))] + #ageFactor = pars.deltaAgeProb[deltaAge(age(man) - age(woman))] # legal dependents (i.e. usually underage persons living at the same house) numChildrenWithWoman = length(dependents(woman)) childrenFactor = 1/exp(pars.bridesChildrenExp * numChildrenWithWoman) - geoFactor * socFactor * ageFactor * childrenFactor * studentFactor + geoFactor * socFactor * ageFactor(age(man), age(woman), pars) * childrenFactor * studentFactor end geoDistance(m, w, pars) = manhattanDistance(getHomeTown(m), getHomeTown(w))/ (pars.mapGridXDimension + pars.mapGridYDimension) -selectMarriage(p, pars) = isMale(p) && isSingle(p) && age(p) > pars.ageOfAdulthood && - careNeedLevel(p) < 4 +selectMarriage(p, pars) = alive(p) && isMale(p) && isSingle(p) && + age(p) > pars.ageOfAdulthood && careNeedLevel(p) < 4 function marriage!(man, time, model, pars) + @assert alive(man) + ageclass = ageClass(man) - manMarriageProb = pars.basicMaleMarriageProb * pars.maleMarriageModifierByDecade[ageclass] + manMarriageProb = ageclass > length(pars.maleMarriageModifierByDecade) ? + 0.0 : pars.basicMaleMarriageProb * pars.maleMarriageModifierByDecade[ageclass] if status(man) != WorkStatus.worker || careNeedLevel(man) > 1 manMarriageProb *= pars.notWorkingMarriageBias end - snc = shareMenNoChildren(model, ageclass) + snc = model.marriageCache.shareMenNoChildren[ageclass+1] den = snc + (1-snc) * pars.manWithChildrenBias prob = manMarriageProb / den * (hasDependents(man) ? pars.manWithChildrenBias : 1) - if rand() >= p_yearly2monthly(prob) + if rand() >= p_yearly2monthly(limit(0.0, prob, 1.0)) return nothing end # get cached list # note: this is getting updated as we go - women = eligibleWomen(model, pars) - - # we store candidates as indices, so that we can efficiently remove married women - candidates = [i for (i,w) in enumerate(women) if (age(man)-10 < age(w) < age(man)+5) && - # exclude siblings as well - !livingTogether(man, w) && !related1stDegree(man, w) ] - - if length(candidates) == 0 + women = model.marriageCache.eligibleWomen + if isempty(women) return nothing end - - weights = [marryWeight(man, women[idx], pars) for idx in candidates] - - cumsum!(weights, weights) + + # keep array across fun calls + weights = model.marriageCache.weights + resize!(weights, length(women)) + sum = 0.0 + for (i,woman) in enumerate(women) + w = marryWeight(man, woman, pars) + weights[i] = sum += w + end + if weights[end] == 0 - selected = rand(1:length(weights)) - else - r = rand() * weights[end] - selected = findfirst(>(r), weights) - @assert selected != nothing + return nothing end - - selectedIdx = candidates[selected] - selectedWoman = women[selectedIdx] + + r = rand() * weights[end] + selected = searchsortedfirst(weights, r) + @assert selected <= length(women) + selectedWoman = women[selected] setAsPartners!(man, selectedWoman) # remove from cached list - remove_unsorted!(women, selectedIdx) + remove_unsorted!(women, selected) joinCouple!(man, selectedWoman, model, pars) diff --git a/src/demography/simulate/relocate.jl b/src/demography/simulate/relocate.jl new file mode 100644 index 0000000..a9f7d38 --- /dev/null +++ b/src/demography/simulate/relocate.jl @@ -0,0 +1,22 @@ +export relocate!, selectRelocate + + +selectRelocate(person, pars) = canLiveAlone(person) && status(person) == WorkStatus.worker && + isSingle(person) && livesInSharedHouse(person) + +function relocate!(person, time, model, pars) + if rand() < pars.moveOutProb + peopleToMove = [person] + for dep in dependents(person) + if livingTogether(person, dep) + push!(peopleToMove, dep) + end + end + + movePeopleToEmptyHouse!(peopleToMove, rand([:near, :far]), model.houses, model.towns) + + return true + end + + false +end diff --git a/src/demography/simulate/socialCareTransition.jl b/src/demography/simulate/socialCareTransition.jl new file mode 100644 index 0000000..a28e0c5 --- /dev/null +++ b/src/demography/simulate/socialCareTransition.jl @@ -0,0 +1,62 @@ +using Distributions +using Utilities + +export selectSocialCareTransition, socialCareTransition! + + +function selectSocialCareTransition(p, pars) + true +end + + +mutable struct SocialCareCache + classBias :: Vector{Float64} +end + +SocialCareCache() = SocialCareCache([]) + + +function socialCarePreCalc!(model, pars) + pc = model.socialCareCache + resize!(pc.classBias, 5) + classes = 0:(length(pars.cumProbClasses)-1) + preCalcRateBias!(classes, pars.careBias, pc.classBias, 1) do c + model.socialCache.socialClassShares[c+1] + end +end + +#=function classSocialCareBias(model, pars, class) + classes = 0:(length(pars.cumProbClasses)-1) + rateBias(classes, pars.careBias, class) do c + model.socialCache.socialClassShares[c+1] + end +end=# + +function socialCareTransition!(person, time, model, pars) + scaling = isFemale(person) ? pars.femaleAgeCareScaling : pars.maleAgeCareScaling + ageCareProb = exp(age(person)/scaling) * pars.personCareProb + + baseProb = pars.baseCareProb + ageCareProb + + class = classRank(person) + if status(person) == WorkStatus.child || status(person) == WorkStatus.student + class = parentClassRank(person) + end + + #baseProb *= classSocialCareBias(model, pars, class) + baseProb *= model.socialCareCache.classBias[class+1] + + if rand() > p_yearly2monthly(limit(0.0, baseProb, 1.0)) + return false + end + + #transitionRate = pars.careTransitionRate * classSocialCareBias(model, pars, class) + transitionRate = pars.careTransitionRate * model.socialCareCache.classBias[class+1] + + careNeed = careNeedLevel(person) + rand(Geometric(1.0-transitionRate)) + 1 + # careLevel is 0-based + careNeedLevel!(person, min(careNeed, numCareLevels(pars)-1)) + true +end + + diff --git a/src/lpm/demography/simulate/socialTransition.jl b/src/demography/simulate/socialTransition.jl similarity index 82% rename from src/lpm/demography/simulate/socialTransition.jl rename to src/demography/simulate/socialTransition.jl index 5e157b5..4a36430 100644 --- a/src/lpm/demography/simulate/socialTransition.jl +++ b/src/demography/simulate/socialTransition.jl @@ -1,8 +1,6 @@ using Distributions: Normal, LogNormal -export socialTransition!, selectSocialTransition - -using XAgents +export socialTransition!, selectSocialTransition function selectSocialTransition(p, pars) @@ -37,27 +35,28 @@ function incomeDist(person, pars) end end -# TODO dummy, replace -# or, possibly remove altogether and calibrate model -# properly instead -socialClassShares(model, class) = 0.2 -function studyClassFactor(person, model, pars) - if classRank(person) == 0 - return socialClassShares(model, 0) > 0.2 ? 1/0.9 : 0.85 - end +mutable struct SocialCache + socialClassShares :: Vector{Float64} +end - if classRank(person) == 1 && socialClassShares(model, 1) > 0.35 - return 1/0.8 - end +SocialCache() = SocialCache([]) - if classRank(person) == 2 && socialClassShares(model, 2) > 0.25 - return 1/0.85 - end +# TODO possibly remove altogether and calibrate model +# properly instead - 1.0 +function socialPreCalc!(model, pars) + pc = model.socialCache + pc.socialClassShares = zeros(5) + + for p in model.pop + pc.socialClassShares[classRank(p)+1] += 1 + end + + pc.socialClassShares ./= length(model.pop) end + doneStudying(person, pars) = classRank(person) >= 4 # TODO @@ -80,11 +79,11 @@ end # probability to start studying instead of working function startStudyProb(person, model, pars) - if father(person) == nothing && mother(person) == nothing + if father(person) == mother(person) == undefinedPerson return 0.0 end - if provider(person) == nothing + if isUndefined(provider(person)) return 0.0 end @@ -102,9 +101,7 @@ function startStudyProb(person, model, pars) (exp(pars.eduWageSensitivity * relCost) + pars.constantIncomeParam) # TODO factor out class - targetEL = father(person) != nothing ? - max(classRank(father(person)), classRank(mother(person))) : - classRank(mother(person)) + targetEL = parentClassRank(person) dE = targetEL - classRank(person) expEdu = exp(pars.eduRankSensitivity * dE) educationEffect = expEdu / (expEdu + pars.constantEduParam) @@ -113,7 +110,7 @@ function startStudyProb(person, model, pars) pStudy = incomeEffect * educationEffect * careEffect - pStudy *= studyClassFactor(person, model, pars) + #pStudy *= studyClassFactor(person, model, pars) return max(0.0, pStudy) end @@ -140,6 +137,8 @@ function startWorking!(person, pars) resetWork!(person, pars) + status!(person, WorkStatus.worker) + dKi = rand(Normal(0, pars.wageVar)) initialIncome!(person, initialIncomeLevel(person, pars) * exp(dKi)) diff --git a/src/generic/XAgents.jl b/src/generic/XAgents.jl deleted file mode 100644 index 6538965..0000000 --- a/src/generic/XAgents.jl +++ /dev/null @@ -1,21 +0,0 @@ -""" -Module for defining a supertype, AbstractAgent for all Agent types - with additional ready-to-use agents to be used in (Multi-)ABMs models -""" -module XAgents - - abstract type AbstractAgent - end - - abstract type AbstractXAgent - end - - getIDCOUNTER() = 0 - - include("../agents/town.jl") - include("../agents/house.jl") - include("../agents/person.jl") - include("../agents/world.jl") - -end # XAgents - diff --git a/src/handleParams.jl b/src/handleParams.jl index 41cf066..5dd4ad5 100644 --- a/src/handleParams.jl +++ b/src/handleParams.jl @@ -6,13 +6,13 @@ using YAML nameOfParType(t) = replace(String(nameof(t)), "Pars" => "") |> Symbol -function saveParametersToFile(simPars::SimulationPars, pars::DemographyPars, fname) +function saveParametersToFile(simPars::SimulationPars, pars::ModelPars, fname) dict = Dict{Symbol, Any}() dict[:Simulation] = parToYaml(simPars) - for f in fieldnames(DemographyPars) - name = nameOfParType(fieldtype(DemographyPars, f)) + for f in fieldnames(ModelPars) + name = nameOfParType(fieldtype(ModelPars, f)) dict[name] = parToYaml(getfield(pars, f)) end @@ -27,8 +27,61 @@ function loadParametersFromFile(fname) simpars = parFromYaml(yaml, SimulationPars, :Simulation) pars = [ parFromYaml(yaml, ft, nameOfParType(ft)) - for ft in fieldtypes(DemographyPars) ] - simpars, DemographyPars(pars...) + for ft in fieldtypes(ModelPars) ] + simpars, ModelPars(pars...) end + +function loadParameters(argv, cmdl...) + arg_settings = ArgParseSettings("run simulation", autofix_names=true) + + @add_arg_table! arg_settings begin + "--par-file", "-p" + help = "parameter file" + default = "" + "--par-out-file", "-P" + help = "file name for parameter output" + default = "parameters.run.yaml" + end + + if ! isempty(cmdl) + add_arg_table!(arg_settings, cmdl...) + end + + # setup command line arguments with docs + + add_arg_group!(arg_settings, "Simulation Parameters") + fieldsAsArgs!(arg_settings, SimulationPars) + + for t in fieldtypes(ModelPars) + groupName = String(nameOfParType(t)) * " Parameters" + add_arg_group!(arg_settings, groupName) + fieldsAsArgs!(arg_settings, t) + end + + # parse command line + args = parse_args(argv, arg_settings, as_symbols=true) + + # read parameters from file if provided or set to default + simpars, pars = loadParametersFromFile(args[:par_file]) + + # override values that were provided on command line + + overrideParsCmdl!(simpars, args) + + @assert typeof(pars) == ModelPars + for f in fieldnames(ModelPars) + overrideParsCmdl!(getfield(pars, f), args) + end + + # pick time-dependent seed if seed == 0 + reseed0!(simpars) + + if args[:par_out_file] != "" + # keep a record of parameters used (including seed!) + saveParametersToFile(simpars, pars, args[:par_out_file]) + end + + simpars, pars, args +end diff --git a/src/lpm/Demography.jl b/src/lpm/Demography.jl deleted file mode 100644 index 3bd59a2..0000000 --- a/src/lpm/Demography.jl +++ /dev/null @@ -1,7 +0,0 @@ -module Demography - - include("./demography/Create.jl") - include("./demography/Initialize.jl") - include("./demography/Simulate.jl") - -end # Demography \ No newline at end of file diff --git a/src/lpm/demography/Initialize.jl b/src/lpm/demography/Initialize.jl deleted file mode 100644 index 9bb7617..0000000 --- a/src/lpm/demography/Initialize.jl +++ /dev/null @@ -1,104 +0,0 @@ -module Initialize - -using Distributions: Normal -using Random: shuffle -using XAgents - -export initializeHousesInTowns, assignCouplesToHouses!, initClass!, initWork! - -"initialize houses in a given set of towns" -function initializeHousesInTowns(towns::Array{Town,1}, pars) - - houses = PersonHouse[] - - for town in towns - if town.density > 0 - - adjustedDensity = town.density * pars.mapDensityModifier - - for hx in 1:pars.townGridDimension - for hy in 1:pars.townGridDimension - - if(rand() < adjustedDensity) - house = PersonHouse(town,(hx,hy)) - push!(houses,house) - end - - end # for hy - end # for hx - - end # if town.density - end # for town - - houses - -end # function initializeHousesInTwons - - -"Randomly assign a population of couples to non-inhebted set of houses" -function assignCouplesToHouses!(population::Array{Person}, houses::Array{PersonHouse}) - women = [ person for person in population if isFemale(person) ] - - randomhouses = shuffle(houses) - - for woman in women - house = pop!(randomhouses) - - moveToHouse!(woman, house) - if !isSingle(woman) - moveToHouse!(partner(woman), house) - end - - for child in dependents(woman) - moveToHouse!(child, house) - end - end # for person - - for person in population - if person.pos == undefinedHouse - @assert isMale(person) - @assert length(randomhouses) >= 1 - moveToHouse!(person, pop!(randomhouses)) - end - end -end # function assignCouplesToHouses - - -function initClass!(person, pars) - p = rand() - class = findfirst(x->p Tables.matrix - deathFemale = CSV.File(deathFFName, header=0) |> Tables.matrix - deathMale = CSV.File(deathMFName, header=0) |> Tables.matrix - - DemographyData(fert,deathFemale,deathMale) -end - -function loadDemographyData(datapars) - dir = datapars.datadir - ukDemoData = loadDemographyData(dir * "/" * datapars.fertFName, - dir * "/" * datapars.deathFFName, - dir * "/" * datapars.deathMFName) -end diff --git a/src/lpm/demography/simulate/birth.jl b/src/lpm/demography/simulate/birth.jl deleted file mode 100644 index c78d7c6..0000000 --- a/src/lpm/demography/simulate/birth.jl +++ /dev/null @@ -1,298 +0,0 @@ -using Utilities - -using XAgents - -export selectBirth, doBirths!, birth! - -function computeBirthProb(rWoman,parameters,data,currstep) - - (curryear,currmonth) = date2yearsmonths(currstep) - currmonth = currmonth + 1 # adjusting 0:11 => 1:12 - - #= - womanClassShares = [] - womanClassShares.append(len([x for x in womenOfReproductiveAge if x.classRank == 0])/float(len(womenOfReproductiveAge))) - womanClassShares.append(len([x for x in womenOfReproductiveAge if x.classRank == 1])/float(len(womenOfReproductiveAge))) - womanClassShares.append(len([x for x in womenOfReproductiveAge if x.classRank == 2])/float(len(womenOfReproductiveAge))) - womanClassShares.append(len([x for x in womenOfReproductiveAge if x.classRank == 3])/float(len(womenOfReproductiveAge))) - womanClassShares.append(len([x for x in womenOfReproductiveAge if x.classRank == 4])/float(len(womenOfReproductiveAge))) - =# - - - if curryear < 1951 - rawRate = parameters.growingPopBirthProb - else - (yearold,tmp) = age2yearsmonths(age(rWoman)) - rawRate = data.fertility[yearold-16,curryear-1950] - end - - #= - a = 0 - for i in range(int(self.p['numberClasses'])): - a += womanClassShares[i]*math.pow(self.p['fertilityBias'], i) - baseRate = rawRate/a - birthProb = baseRate*math.pow(self.p['fertilityBias'], womanRank) - =# - - # The above formula with one single socio-economic class translates to: - - birthProb = rawRate * parameters.fertilityBias - return birthProb -end # computeBirthProb - - -function effectsOfMaternity!(woman, pars) - startMaternity!(woman) - - workingHours!(woman, 0) - income!(woman, 0) - potentialIncome!(woman, 0) - availableWorkingHours!(woman, 0) - # commented in sim.py: - # woman.weeklyTime = [[0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12, [0]*12+[1]*12] - # sets all weeklyTime slots to 1 - # TODO copied from the python code, but does it make sense? - setFullWeeklyTime!(woman) - #= TODO - woman.maxWeeklySupplies = [0, 0, 0, 0] - woman.residualDailySupplies = [0]*7 - woman.residualWeeklySupplies = [x for x in woman.maxWeeklySupplies] - =# - - # TODO not necessarily true in many cases - if provider(woman) == nothing - setAsProviderProvidee!(partner(woman), woman) - end - - nothing -end - - -function birth!(woman, currstep, data, parameters) - - # womanClassRank = woman.classRank - # if woman.status == 'student': - # womanClassRank = woman.parentsClassRank - - birthProb = computeBirthProb(woman, parameters, data, currstep) - - assumption() do - @assert isFemale(woman) - @assert ageYoungestAliveChild(woman) > 1 - @assert !isSingle(woman) - @assert age(woman) >= parameters.minPregnancyAge - @assert age(woman) <= parameters.maxPregnancyAge - @assert birthProb >= 0 - end - - #= - The following code is commented in the python code: - #baseRate = self.baseRate(self.socialClassShares, self.p['fertilityBias'], rawRate) - #fertilityCorrector = (self.socialClassShares[woman.classRank] - self.p['initialClassShares'][woman.classRank])/self.p['initialClassShares'][woman.classRank] - #baseRate *= 1/math.exp(self.p['fertilityCorrector']*fertilityCorrector) - #birthProb = baseRate*math.pow(self.p['fertilityBias'], woman.classRank) - =# - - if rand() < p_yearly2monthly(birthProb) - - # parentsClassRank = max([woman.classRank, woman.partner.classRank]) - # baby = Person(woman, woman.partner, self.year, 0, 'random', woman.house, woman.sec, -1, - # parentsClassRank, 0, 0, 0, 0, 0, 0, 'child', False, 0, month) - - baby = Person(pos=woman.pos, - father=partner(woman),mother=woman, - gender=rand([male,female])) - - # this goes first, so that we know material circumstances - effectsOfMaternity!(woman, parameters) - - setAsGuardianDependent!(woman, baby) - if !isSingle(woman) # currently not an option - setAsGuardianDependent!(partner(woman), baby) - end - setAsProviderProvidee!(woman, baby) - - return baby - end # if rand() - - nothing -end - -function verboseBirthCounting(people,parameters) - - allFemales = [ woman for woman in people if isFemale(woman) ] - adultWomen = [ aWoman for aWoman in allFemales if - age(aWoman) >= parameters.minPregnancyAge ] - notFertiledWomen = [ nfWoman for nfWoman in adultWomen if - age(nfWoman) > parameters.maxPregnancyAge ] - womenOfReproductiveAge = [ rWoman for rWoman in adultWomen if - age(rWoman) <= parameters.maxPregnancyAge ] - marriedWomenOfReproductiveAge = - [ rmWoman for rmWoman in womenOfReproductiveAge if - !isSingle(rmWoman) ] - womenWithRecentChild = [ rcWoman for rcWoman in adultWomen if - ageYoungestAliveChild(rcWoman) <= 1 ] - reproductiveWomen = [ rWoman for rWoman in marriedWomenOfReproductiveAge if - ageYoungestAliveChild(rWoman) > 1 ] - womenOfReproductiveAgeButNotMarried = - [ rnmWoman for rnmWoman in womenOfReproductiveAge if - isSingle(rnmWoman) ] - - # for person in self.pop.livingPeople: - # - # if person.sex == 'female' and person.age >= self.p['minPregnancyAge']: - # adultLadies += 1 - # if person.partner != None: - # marriedLadies += 1 - # marriedPercentage = float(marriedLadies)/float(adultLadies) - - numMarriedRepLadies = length(womenOfReproductiveAge) - - length(womenOfReproductiveAgeButNotMarried) - repMarriedPercentage = numMarriedRepLadies / length(adultWomen) - womenWithRecentChildPercentage = length(womenWithRecentChild) / numMarriedRepLadies - - println("# allFemales : $(length(allFemales))") - println("# adult women : $(length(adultWomen))") - println("# NotFertile : $(length(notFertiledWomen))") - println("# fertile women : $(length(womenOfReproductiveAge))") - println("# non-married fertile women : $(length(womenOfReproductiveAgeButNotMarried))") - println("# of women with recent child: $(length(womenWithRecentChild))") - println("married reproductive percentage : $repMarriedPercentage") - println(" out of which had a recent child : $womenWithRecentChildPercentage ") - - nothing -end -""" -Accept a population and evaluates the birth rate upon computing -- the population of married fertile women according to -fixed parameters (minPregnenacyAge, maxPregnenacyAge) and -- the birth probability data (fertility bias and growth rates) - -Class rankes and shares are temporarily ignored. -""" - -selectBirth(woman, parameters) = isFemale(woman) && - !isSingle(woman) && - age(woman) >= parameters.minPregnancyAge && - age(woman) <= parameters.maxPregnancyAge && - ageYoungestAliveChild(woman) > 1 - - -function doBirths!(people, currstep, data, parameters) - - assumption() do - for person in people - @assert alive(person) - end - end - - babies = Person[] - # numBirths = 0 # instead of [0, 0, 0, 0, 0] - - # TODO The following could be collapsed into one loop / not sure if it is more efficient - # there is also a potential to save alot of re-computation in each iteration by - # storing the intermediate results and modifying the computation. - # However, it could be also the case that Julia compiler does something efficient any way? - - reproductiveWomen = [ woman for woman in people if selectBirth(woman, parameters) ] - - # TODO @assumption - assumption() do - allFemales = [ woman for woman in people if isFemale(woman) ] - adultWomen = [ aWomen for aWomen in allFemales if - age(aWomen) >= parameters.minPregnancyAge ] - nonadultFemale = setdiff(Set(allFemales),Set(adultWomen)) - for woman in nonadultFemale - @assert(isSingle(woman)) - @assert !hasChildren(woman) - end - - for woman in allFemales - if woman ∉ reproductiveWomen - @assert isSingle(woman) || - age(woman) < parameters.minPregnancyAge || - age(woman) > parameters.maxPregnancyAge || - ageYoungestAliveChild(woman) <= 1 - end - end - end - - delayedVerbose() do - (curryear,currmonth) = date2yearsmonths(currstep) - currmonth += 1 # adjusting 0:11 => 1:12 - - # TODO this generic print msg to be placed in a top function - println("In iteration $curryear , month $currmonth :") - verboseBirthCounting(people,parameters) - end # verbose - - - #= - adultLadies_1 = [x for x in adultWomen if x.classRank == 0] - marriedLadies_1 = len([x for x in adultLadies_1 if x.partner != None]) - if len(adultLadies_1) > 0: - marriedPercentage.append(marriedLadies_1/float(len(adultLadies_1))) - else: - marriedPercentage.append(0) - adultLadies_2 = [x for x in adultWomen if x.classRank == 1] - marriedLadies_2 = len([x for x in adultLadies_2 if x.partner != None]) - if len(adultLadies_2) > 0: - marriedPercentage.append(marriedLadies_2/float(len(adultLadies_2))) - else: - marriedPercentage.append(0) - adultLadies_3 = [x for x in adultWomen if x.classRank == 2] - marriedLadies_3 = len([x for x in adultLadies_3 if x.partner != None]) - if len(adultLadies_3) > 0: - marriedPercentage.append(marriedLadies_3/float(len(adultLadies_3))) - else: - marriedPercentage.append(0) - adultLadies_4 = [x for x in adultWomen if x.classRank == 3] - marriedLadies_4 = len([x for x in adultLadies_4 if x.partner != None]) - if len(adultLadies_4) > 0: - marriedPercentage.append(marriedLadies_4/float(len(adultLadies_4))) - else: - marriedPercentage.append(0) - adultLadies_5 = [x for x in adultWomen if x.classRank == 4] - marriedLadies_5 = len([x for x in adultLadies_5 if x.partner != None]) - if len(adultLadies_5) > 0: - marriedPercentage.append(marriedLadies_5/float(len(adultLadies_5))) - else: - marriedPercentage.append(0) -=# - - for woman in reproductiveWomen - - baby = birth!(woman, currstep, data, parameters) - if baby != nothing - push!(babies,baby) - end - - end # for woman - - delayedVerbose() do - println("number of births : $length(babies)") - end - - # any reason for that? -# return (newbabies=babies) - - babies -end # function doBirths! - -"This function is supposed to implement the suggested model, TODO" -function doBirthsOpt() end - -# the following accessory functions to be moved to an internal module -population(model) = model.pop -data(model) = model -alivePeople(model) = Iterators.filter(a->alive(a), population(model)) -birthPars(pars) = pars.birthpars - -# Generic API for doDeaths! -doBirths!(model,time,parameters) = - doBirths!(alivePeople(model),time,data(model),birthPars(parameters)) - - - - - diff --git a/src/lpm/demography/simulate/death.jl b/src/lpm/demography/simulate/death.jl deleted file mode 100644 index ec15f7b..0000000 --- a/src/lpm/demography/simulate/death.jl +++ /dev/null @@ -1,181 +0,0 @@ - -using Utilities: age2yearsmonths, date2yearsmonths - -using XAgents: Person, isMale, isFemale, alive -using XAgents: age - -export doDeaths!, setDead! - -function deathProbability(baseRate,person,parameters) - #= - Not realized yet / to be realized in another module? - classRank = person.classRank - if person.status == 'child' or person.status == 'student': - classRank = person.parentsClassRank - =# - - if isMale(person) - mortalityBias = parameters.maleMortalityBias - else - mortalityBias = parameters.femaleMortalityBias - end - - #= - To be integrated in class modules - a = 0 - for i in range(int(self.p['numberClasses'])): - a += self.socialClassShares[i]*math.pow(mortalityBias, i) - =# - - #= - if a > 0: - lowClassRate = baseRate/a - classRate = lowClassRate*math.pow(mortalityBias, classRank) - deathProb = classRate - - b = 0 - for i in range(int(self.p['numCareLevels'])): - b += self.careNeedShares[classRank][i]*math.pow(self.p['careNeedBias'], (self.p['numCareLevels']-1) - i) - - if b > 0: - higherNeedRate = classRate/b - deathProb = higherNeedRate*math.pow(self.p['careNeedBias'], (self.p['numCareLevels']-1) - person.careNeedLevel) # deathProb - =# - - # assuming it is just one class and without care need, - # the above code translates to: - - deathProb = baseRate * mortalityBias - - ##### Temporarily by-passing the effect of Unmet Care Need ############# - - # The following code is already commented in the python code - # a = 0 - # for x in classPop: - # a += math.pow(self.p['unmetCareNeedBias'], 1-x.averageShareUnmetNeed) - # higherUnmetNeed = (classRate*len(classPop))/a - # deathProb = higherUnmetNeed*math.pow(self.p['unmetCareNeedBias'], 1-shareUnmetNeed) - deathProb -end # function deathProb - - -function setDead!(person) - person.info.alive = false - resetHouse!(person) - if !isSingle(person) - resolvePartnership!(partner(person),person) - end - - # dead persons are no longer dependents - setAsIndependent!(person) - - # dead persons no longer have to be provided for - setAsSelfproviding!(person) - - for p in providees(person) - provider!(p, nothing) - # TODO update provision/work status - end - empty!(providees(person)) - - # dependents are being taken care of by assignGuardian! - nothing -end - - -# currently leaves dead agents in population -function death!(person, currstep, data, parameters) - - (curryear,currmonth) = date2yearsmonths(currstep) - currmonth += 1 # adjusting 0:11 => 1:12 - - agep = age(person) - - assumption() do - @assert alive(person) - @assert isMale(person) || isFemale(person) # Assumption - @assert typeof(agep) == Rational{Int} - end - - if curryear >= 1950 - - agep = agep > 109 ? Rational(109) : agep - ageindex = trunc(Int,agep) - rawRate = isMale(person) ? data.deathMale[ageindex+1,curryear-1950+1] : - data.deathFemale[ageindex+1,curryear-1950+1] - - # lifeExpectancy = max(90 - agep, 3 // 1) # ??? This is a direct translation - - else # curryear < 1950 / made-up probabilities - - babyDieProb = agep < 1 ? parameters.babyDieProb : 0.0 # does not play any role in the code - ageDieProb = isMale(person) ? - exp(agep / parameters.maleAgeScaling) * parameters.maleAgeDieProb : - exp(agep / parameters.femaleAgeScaling) * parameters.femaleAgeDieProb - rawRate = parameters.baseDieProb + babyDieProb + ageDieProb - - # lifeExpectancy = max(90 - agep, 5 // 1) # ??? Does not currently play any role - - end # currYear < 1950 - - #= - Not realized yet - classPop = [x for x in self.pop.livingPeople - if x.careNeedLevel == person.careNeedLevel] - Classes to be considered in a different module - =# - - deathProb = min(1.0, deathProbability(rawRate,person,parameters)) - - #= - The following is uncommented code in the original code < 1950 - #### Temporarily by-passing the effect of unmet care need ###### - # dieProb = self.deathProb_UCN(rawRate, person.parentsClassRank, person.careNeedLevel, person.averageShareUnmetNeed, classPop) - =# - - if rand() < p_yearly2monthly(deathProb) - delayedVerbose() do - y, m = age2yearsmonths(agep) - println("person $(person.id) died year $(curryear) with age of $y") - end - setDead!(person) - return true - # person.deadYear = self.year - # deaths[person.classRank] += 1 - end # rand - - false -end - - - -"evaluate death events in a population" -function doDeaths!(people, currstep, data, parameters) - - deads = Person[] - - for person in people - if death!(person, currstep, data, parameters) - push!(deads,person) - end - end # for livingPeople - - delayedVerbose() do - count = length([person for person in people if alive(person)] ) - numDeaths = length(deads) - println("# living people : $(count+numDeaths), # people died in curr iteration : $(numDeaths)") - end - - deads -end # function doDeaths! - -# Atiyah : Alternative: - -# the following accessory functions to be moved to an internal module -population(model) = model.pop -alivePeople(model) = Iterators.filter(a->alive(a), population(model)) -populationPars(pars) = pars.poppars - -# Generic API for doDeaths! -doDeaths!(model,time,parameters) = - doDeaths!(alivePeople(model),time,model,populationPars(parameters)) \ No newline at end of file diff --git a/src/malpm/Demography.jl b/src/malpm/Demography.jl deleted file mode 100644 index a72f295..0000000 --- a/src/malpm/Demography.jl +++ /dev/null @@ -1,52 +0,0 @@ -module Demography - -using XAgents: Town, PersonHouse, Person -using MultiAgents: AbstractMABM, ABM -# using LPM.Demography.Create: createTowns - -export MAModel - -import MultiAgents.Util: AbstractExample -import MultiAgents: allagents -export allagents, population -export DemographyExample, LPMUKDemography, LPMUKDemographyOpt - -### Example Names -"Super type for all demographic models" -abstract type DemographyExample <: AbstractExample end - -"This corresponds to direct translation of the python model" -struct LPMUKDemography <: DemographyExample end - -"This is an attemp for improved algorthimic translation" -struct LPMUKDemographyOpt <: DemographyExample end - - -mutable struct MAModel <: AbstractMABM - towns :: ABM{Town} - houses :: ABM{PersonHouse} - pop :: ABM{Person} - - function MAModel(model,pars,data) - ukTowns = ABM{Town}(model.towns,parameters = pars.mappars) - ukHouses = ABM{PersonHouse}(model.houses) - parameters = (poppars = pars.poppars, birthpars = pars.birthpars, - divorcepars = pars.divorcepars, workpars = pars.workpars) - ukPopulation = ABM{Person}(model.pop,parameters=pars,data=data) - new(ukTowns,ukHouses,ukPopulation) - end - -end - -allagents(model::MAModel) = allagents(model.pop) -population(model::MAModel) = allagents(model.pop) -houses(model::MAModel) = allagents(model.houses) -towns(model::MAModel) = allagents(model.towns) - - -include("./demography/Population.jl") -include("./demography/Simulate.jl") -include("./demography/SimSetup.jl") - - -end # Demography \ No newline at end of file diff --git a/src/malpm/demography/Population.jl b/src/malpm/demography/Population.jl deleted file mode 100644 index ebae63f..0000000 --- a/src/malpm/demography/Population.jl +++ /dev/null @@ -1,74 +0,0 @@ -""" -Population module providing help utilities for realizing a population as an ABM -""" - -module Population - -using MultiAgents: ABM, AbstractABMSimulation, AbstractMABM -using MultiAgents: allagents, dt, kill_agent! -using XAgents: Person -using XAgents: alive, agestepAlive! -# using MALPM.Demography: population - -import XAgents: agestep! - -export population_step!, agestepAlivePerson!, removeDead! - - -"Step function for the population" -function population_step!(population::ABM{PersonType}, - sim::AbstractABMSimulation, - example) where {PersonType} - for person in allagents(population) - if alive(person) - agestep!(person,dt(sim)) - end - end -end - -population_step!(model::AbstractMABM, - sim::AbstractABMSimulation, - example) = - population_step!(model.pop,sim,example) - -"remove dead persons" -function removeDead!(person::PersonType, population::ABM{PersonType}) where {PersonType} - if !alive(person) - kill_agent!(person, population) - end - nothing -end - -function removeDead!(population::ABM{PersonType}, - simulation::AbstractABMSimulation, - example) where {PersonType} - people = reverse(allagents(population)) - for person in people - alive(person) ? nothing : kill_agent!(person,population) - end - nothing -end - -"increment age with the simulation step size" -agestep!(person::Person,population::ABM{Person}, - sim::AbstractABMSimulation, - example) where {PersonType} = agestep!(person,dt(sim)) - - -agestep!(person::Person,model::AbstractMABM,sim,example) = - agestep!(person,model.pop,sim,example) - -"increment age with the simulation step size" -agestepAlivePerson!(person::PersonType,population::ABM{PersonType}, - sim::AbstractABMSimulation, - example) where {PersonType} = - agestepAlive!(person, dt(sim)) - - -agestepAlivePerson!(person,model::AbstractMABM,sim,example) = - agestepAlivePerson!(person,model.pop,sim,example) - -end # Population - - - diff --git a/src/malpm/demography/SimSetup.jl b/src/malpm/demography/SimSetup.jl deleted file mode 100644 index 2ff0e04..0000000 --- a/src/malpm/demography/SimSetup.jl +++ /dev/null @@ -1,61 +0,0 @@ -module SimSetup - - -using MALPM.Demography.Population: removeDead!, - agestepAlivePerson!, agestep!, population_step! - -using MALPM.Demography: DemographyExample, - LPMUKDemography, LPMUKDemographyOpt -using MALPM.Demography.Simulate: doDeaths!, doBirths!, doDivorces! - -using MultiAgents: AbstractABMSimulation -using MultiAgents: attach_pre_model_step!, attach_post_model_step!, - attach_agent_step! -using Utilities: setVerbose!, unsetVerbose!, setDelay!, - checkAssumptions!, ignoreAssumptions! -import MultiAgents: setup!, verbose -export setup! - -""" -set simulation paramters @return dictionary of symbols to values - -All information needed by the generic Simulations.run! function -is provided here - -@return dictionary of required simulation parameters -""" - - -function setupCommon!(sim::AbstractABMSimulation) - - verbose(sim) ? setVerbose!() : unsetVerbose!() - setDelay!(sim.parameters.sleeptime) - sim.parameters.checkassumption ? checkAssumptions!() : - ignoreAssumptions!() - - attach_post_model_step!(sim,doDeaths!) - attach_post_model_step!(sim,doBirths!) - attach_post_model_step!(sim,doDivorces!) - nothing -end - -"set up simulation functions where dead people are removed" -function setup!(sim::AbstractABMSimulation, example::LPMUKDemography) - # attach_pre_model_step!(sim,population_step!) - attach_agent_step!(sim,agestep!) - setupCommon!(sim) - - nothing -end - - -function setup!(sim::AbstractABMSimulation,example::LPMUKDemographyOpt) - - attach_agent_step!(sim,agestepAlivePerson!) - setupCommon!(sim) - - nothing -end - - -end # SimSetup \ No newline at end of file diff --git a/src/malpm/demography/Simulate.jl b/src/malpm/demography/Simulate.jl deleted file mode 100644 index 28e5c95..0000000 --- a/src/malpm/demography/Simulate.jl +++ /dev/null @@ -1,88 +0,0 @@ -""" - Main simulation functions for the demographic aspect of LPM. -""" - -module Simulate - -# using MultiAgents.Util: getproperty - -using XAgents: Person, isFemale, alive, age - -using MultiAgents: ABM, AbstractMABM, AbstractABMSimulation -using MultiAgents: allagents, add_agent!, currstep, verbose -using MALPM.Demography.Population: removeDead! -using MALPM.Demography: DemographyExample, LPMUKDemography, LPMUKDemographyOpt, - houses, towns -using LPM -import LPM.Demography.Simulate: doDeaths!, doBirths!, doDivorces! -# export doDeaths!,doBirths! - -alivePeople(population::ABM{Person},::LPMUKDemography) = allagents(population) - -alivePeople(population::ABM{Person},::LPMUKDemographyOpt) = - # Iterators.filter(person->alive(person),allagents(population)) - [ person for person in allagents(population) if alive(person) ] - -function removeDeads!(deadpeople,population,::LPMUKDemography) - for deadperson in deadpeople - removeDead!(deadperson,population) - end - - nothing -end - -removeDeads!(deadpeople,population,::LPMUKDemographyOpt) = nothing - -function doDeaths!(model::AbstractMABM, sim::AbstractABMSimulation, example::DemographyExample) # argument simulation or simulation properties ? - - population = model.pop - - (deadpeople) = LPM.Demography.Simulate.doDeaths!( - alivePeople(population,example), - currstep(sim), - population.data, - population.parameters.poppars) - - removeDeads!(deadpeople,population,example) - nothing -end # function doDeaths! - - -function doBirths!(model::AbstractMABM, sim::AbstractABMSimulation, example::DemographyExample) - - population = model.pop - - newbabies = LPM.Demography.Simulate.doBirths!( - alivePeople(population,example), - currstep(sim), - population.data, - population.parameters.birthpars) - - # false ? population.variables[:numBirths] += length(newbabies) : nothing # Temporarily this way till realized - - for baby in newbabies - add_agent!(population,baby) - end - - nothing -end - - -function doDivorces!(model::AbstractMABM, sim::AbstractABMSimulation, example::DemographyExample) - - population = model.pop - - divorced = LPM.Demography.Simulate.doDivorces!( - allagents(population), - currstep(sim), - houses(model), - towns(model), - population.parameters.divorcepars) - - nothing -end - - - - -end # Simulate \ No newline at end of file diff --git a/src/multiagents/XAgents.jl b/src/multiagents/XAgents.jl deleted file mode 100644 index b0130d9..0000000 --- a/src/multiagents/XAgents.jl +++ /dev/null @@ -1,15 +0,0 @@ -""" -Module for defining a supertype, AbstractAgent for all Agent types - with additional ready-to-use agents to be used in (Multi-)ABMs models -""" -module XAgents - - using MultiAgents: AbstractAgent, AbstractXAgent, getIDCOUNTER - - include("../agents/town.jl") - include("../agents/house.jl") - include("../agents/person.jl") - include("../agents/world.jl") - -end # XAgents -