Creating Scripts
=================

In this section, we explore ``Cytocalc`` by creating two simple scripts that uses some functions in built in ``Cytocalc``. The :doc:`code
reference</cytocalc>` lists all functions available for usage and some example scripts can
be found in the ``scripts`` directory.

General idea
------------

#. Import the specific class of ``Cytocalc`` you want to use and additional python modules.
#. Define an object of that specific class
#. Read over data stored in single or multiple files and add it to the object
#. Call the specific analysis function you want to use
#. Plot or print result 

Example 1: MSD of beads
-----------------------

**Simulation details:** A cytosim simulation with *N* beads in a filament network.
The report is then generated using ``report bead:position beads.txt``. 

#. We start by importing the :class:`.CSMParser` for data analysis, and ``matplotlib`` for ploting. 

        .. code-block:: py

                import matplotlib.pyplot as plt
                from cytocalc.csmparser import CSMParser

#. We want to read the report file with the name ``<path-to-beads.txt>``.

        .. code-block:: py

                filename = '<path-to-beads.txt>'

#. We then use the :meth:`.parse_simFile` function of the :class:`.CSMParser` to read the report file. We use the :meth:`.parse_simFile` function  because the data from every frame is contained in a single text file.

        .. code-block:: py

                parser = CSMParser()
                sim = parser.parse_simFile(filename)

#. We can now use the :meth:`.compute_time_avg_msd` function to compute the Mean square displacement of the beads in the system.

        .. code-block:: py

                msd = sim.compute_time_avg_msd()
               
#. We use ``matplotlib`` to plot the outputs of the :meth:`.compute_time_avg_msd` function, which shows time vs mean square displacement of the beads. 

        .. code-block:: py

                plt.plot(msd.keys(),msd.values())
                plt.show()                

Example 2: Contraction of network
----------------------------------

**Simulation details:** Following the studies on effects of crosslinkers and motors on the contraction of cytoskeleton network [1]_, here we have a set of cytosim simulations of filament networks with motors and crosslinkers. A parameter scan over different number of motors is saved in the folders named ``motor_count####``, and the position of each filament nodes for each frame of the simulation are reported using ``reportf``, creating a directory structure as follows:

.. code-block:: sh

        .
        ├── config.cym
        ├── motor_count1000
        │   ├── report0001.txt
        │   ├── ...
        │   └── report0040.txt
        ├── motor_count2000
        │   ├── report0001.txt
        │   ├── ...
        │   └── report0040.txt
        ├── motor_count3000
        │   ├── report0001.txt
        │   ...
        ...



#. We import the ``glob`` module for finding certain files in a directory. From our package, we also import the CSMSimulation class.

        .. code-block:: py

                import glob
                from cytocalc.csmsimulation import CSMSimulation

#. We get the name of the directory, ``dir`` containing a few ``report####.txt`` files corresponding to the frames of a simulation. The directory name is provided as the first argument when calling the python script. We then find the `.txt` files in this directory using the ``glob`` module.

        .. code-block:: py

                dir = sys.argv[1]
                frame_files = glob.glob(dir+'*.txt')

#. We read each of the files corresponding to individual frames using the :meth:`.parse_frameFile` function and add it to a :class:`.CSMSimulation`.
        
        .. code-block:: py

                sim = CSMSimulation()
                for frame_file in frame_files:
                        frame = CSMParser.parse_frameFile(frame_file)
                        sim.add_frame(frame)

#. We compute the radius of network in each frame using the :meth:`.compute_radius` function.

        .. code-block:: py

                times = [frame.csm.time for frame in sim.frames]
                radii = [frame.compute_radius() for frame in sim.frames]
        

#. We can plot  radius of network w.r.t time .

        .. code-block:: py
                
                plt.plot(times, radii)   
        .. image:: /plots/radii_vs_time_750.png
                :align: center
                :width: 400

#. We can also compute the mean contraction rate of the network using the :meth:`.contraction_rate` function.

        .. code-block:: py

                print(sim.contraction_rate())

        This will return the contraction rate in Cytosim units (μm/s). 

        .. code-block:: sh

                $ ./contrate.py ../data/2022_01_10_22_17_56/
                -0.07988585749886662

#. Looping over all directories ``motor_count####``, we can also observe how contraction rate varies with different parameters, which is motor count in this particular case, along with the thoretical contraction rate [1]_ plotted using :func:`theoretical_cont_rate`.

        .. image:: /plots/crate_vs_motors.png
                :align: center
                :width: 400

.. [1] Belmonte et. al. 2017, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5615920/